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Abstract 


In almost all numerical simulations matrix inversion is the most recursive and 
computationally intensive step Therefore, the choice of matrix inveitei , or equiv- 
alently the linear system solver is a ciitical step towards reducing the compu- 
tational time of the o^elall simulation Conjugate Giadient method is an ex- 
tremely effective technique for sohing linear systems when the coefficient matrix 
is symmetiic and positive definite Generalization and efficient implementation 
of conjugate giadient method foi the linear systems with unsymmetric coefficient 
matiices is still an active area of leseaich One such generalization, based on pre- 
conditioning with incomplete lowei-upper (ILU) decomposition of the coefficient 
matrix, has been proposed by Kershaw [1978] Following this generalization, She- 
orey [2001] has developed a Preconditioned Conjugate Gradient (PCG) algorithm 
that IS believed to w'oik well for the case of any nonsmgular coefficient matrix 
The aim of the present leseaich is to assess the above mentioned PCG algo- 
rithm for matrix inversion in various test cases including the complex problem 
of Czochralski crystal growth simulation Czochialski technique is quite popular 
foi growing single crystals that are of great importance for modern electronic and 
electro-optical applications and devices Performance of these devices is deleteri- 
ouslv influenced by the presence of solute inhomogeneities m the crystals used In 
order to investigate the solute segregation m these crystals, the coupled problem 
of heat transfer, hydiodynamics and solute consenation must be simultaneously 
solved This research presents a macroscopic model for solute segregation 

CPU time analysis for the test case of unsteady-nonlinear heat conduction 
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shows that PCG is e\tiemely effective for matiix inversion as compaied to Gaus- 
sian elimination and Gauss-Seidel iterations PCG becomes an e^en moic at- 
tractive choice on finci grids Fuithei PCG is found to be well suited foi matrix 
m/eision in advection-diffusion problems The results obtained loi the fully devel- 
oped Nusselt number in channel flow matches closely with the analytical solution 
lepoited m the hteiatuie A compaiison of the eigenvalue spectra of the oiigmal 
and preconditioned matrices show that pieconditiomng with incomplete lowei- 
upper (ILU) decomposition is very effectne m lowering the condition number 
towards unity Moreover, this preconditioning technique is general m the sense 
that it can be constructed and applied to any nonsingulai coefficient matrix 
For growth of a Nd YAG crystal from its melt, the PCG algorithm has been 
used to determine the quasi-steady flow and temperatuie fields In addition, the 
distribution of Nd particles m YAG has been determined by a truly unsteady 
calculation Results obtained in this part of the study show the following 

In the simulations of Czochralski ci;>stal growth, the pull velocity is varied 
with time to grow constant diameter crystals The results obtained foi pull- 
velocity variation show that it is not possible to maintain favorable growth con- 
ditions over a longer growth period This is because the pull velocity continu- 
ously decreases with time After a small inciease in crystal length, melting starts 
Therefore, to grow longer crystals with constant diameter, other parameters such 
as the crucible temperature and the enclosure temperature need to be varied, 
while keeping pull velocity positive and a constant Redistribution of the seg- 
regated solute is found to correlate with the complex flow patterns iii the melt 
Profiles of the radial solute distribution reveal that crystal rotation improves ra- 
dial uniformity of solute concentration m the gioun crystals, while the uniformity 
along the crystal length is unaffected 
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Chapter 1 
Introduction 


In all numerical simulations, the lesult of the disci etization process is a system 
of algebraic equations, which can be linear or non-linear according to the nature 
of the partial differential equation(s) that govern the physical process of interest 
In the non-lineai case, the discretized equations must be sohed by an iterative 
technique that involve guessing a solution, linearizing the equations about that 
solution, and improving the solution until a comerged result is obtained More- 
over, if the mathematical problem is unsteady, these nonlinearity iterations need 
to be carried out at every time step Thus in numerical simulations, matrix in- 
version is the most recursive and computationally intensive step and therefore, 
the choice of matrix inverter, or equivalently the linear system solver is a critical 
step towards reducing the computational time of the overall simulation 

Matrix inverters can be classified as direct and iterative Direct methods are 
usually a variation of the Gaussian elimination technique, making use of forward 
elimination and backward substitution It is equivalent to the decomposition of 
the coefficient matrix into lower and upper triangular factors Diiect inverters in- 
volve larger number of arithmetic operations per inversion and hence carry large 
round-off errors Moreover, direct sobers when applied to a sparse system, in- 
troduce a large number of additional nonzero entries into the coefficient matrix 
These ^'fill-in" values destroy the sparse structure of the coefficient matrix and 
thus increases the storage lequiiements significantly Howe^el direct methods 
obtain the solution in finite number of predetermined steps and are generally 
stable On the other hand, the term ‘ iterative methods" refers to a wide range of 
techniques that use successive appioximations to obtain more accurate solutions 
to a linear system at each step The computational effort per inversion of an 
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iterative technique depends on its convergence rate The rate of convergence of 
these methods depends gieatly on the eigenvalue spectrum of the coefficient ma- 
trix Hence, iterative methods usually involve a second matrix that transforms 
the coefficient matrix into one w'th a more favorable spectrum This technique 
IS called '^preconditioning" and the transformation matrix is called a “precondi- 
tioned' A good pieconditionei improves the convergence rate of an rteratrve 
method, sufficiently to overcome the extra cost of constructing and applying the 
pieconditionei When using an iterative method, the coefficient matrix is in- 
\olved (directly or indirectly) only m terms of matrix by vector products Thus 
there is no need for storing the zero elements of the coefficient matrix, and there- 
fore, it is possible to implement the algorithms that are extremely cost-effective 
with respect to computing time as well as storage In addition of the matrix by 
vector product, many iterative solvers are based on the vector operations like the 
calculation of a inner product Thus except for some possible problems related to 
the global communication needed for the computation of inner products, such it- 
erative methods are well suited for parallel and vectorized implementations But 
then iterative methods are only conditionally stable To combine the benefits of 
both the direct and iterative approaches, hybrid methods have been developed 
The “Preconditioned Conjugate Gradient (PCG)" is one such technique that has 
been assessed m the present study 

1.1 Literature Review on Matrix Inversion 

Lanczos [1952] and Hestenes and Stiefel [1952] discovered the method of Con- 
jugate Gradients (CG) for solving linear systems with symmetric and positive 
definite coefficient matrix A This method was initially regarded as an exact (di- 
rect) solution method, guaranteed to converge in no more than n steps, where n 
IS the dimension of the problem Moie precisely, m exact arithmetic the method 
is guaranteed to compute the solution to Ax = b m exactly m steps, wheie m is 
the degree of minimal polynomial ot A, i e , the number of distinct eigenvalues of 
A (since A is assumed to be symmetiic) Based on numerical experiments it was 
observed that in presence of rounding errors, a few more than n iterations were 
often necessary to obtain an accurate solution, paiticulaily on ill-conditioned 
systems On the other hand, Hestenes and Stiefel [1952] observed that on well 
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conditioned systems, satisfactoi}^ appioximations could be obtained in far fewer 
than n iterations, making the Conjugate Giadient method an attractive alter- 
native to Gaussian Elimination In spite of this, for a number of leasons the 
Conjugate Giadient method saw lelativelj^ little use for almost 20 yeais, during 
which Gaussian Elimination (foi dense matiices) and SOR and Chebyschev iter- 
ation (for sparse matiices) weie usually the preferred solution methods Great 
details of these developments can be found in the interesting paper on the history 
of conjugate gradient method by Golub and O’Leary [1989] 

This situation finally changed aiound 1970 with the publication of a paper 
of Reid Reid [1971] shown that for large sparse matrices that aie reasonably 
well conditioned, the Conjugate Gradient method is a very powerful scheme that 
yields good approximate solutions m far fewer than n steps Reid’s paper set the 
stage for the rehabilitation of the Conjugate Gradient method and stimulated 
research m ^arlous directions On one hand, extensions to linear systems with 
indefinite and/or nons\mmetiic matrices were sought, on the other, techniques 
to improve the conditioning of linear systems were developed m order to improve 
the rate of convergence of the method 

Early ideas of preconditioning the Conjugate Gradient method can be found 
m Engeh et al [1959], and occasionally in several papers published during 1960s 
(Golub and O’Leary [1989]) A detailed study of SSOR preconditioning as a 
means of accelerating the convergence was carried out in by Axelsson [1972] 

A measure breakthrough took place around the mid-1970s, with the intro- 
duction by Meijerink and van der Vorst of the incomplete Cholesky-Conjugate 
Gradient (ICCG) algorithm Incomplete factoiwation methods were introduced 
for the first time by Buleev m the then-Soviet Union m the late 1950s, and inde- 
pendently by Varga m 1970s However, Meijeiink and van der Vorst [1977] deserve 
credit for recognizing the potential of incomplete factorizations as preconditioners 
for the Conjugate Gradient method A.nothei paper that did much to popularize 
this approach is by Kershaw [1978] Since then, a number of improvements and 
extensions have been made, including level-of-fills and drop tolerances-based in- 
complete factorizations, generalizations to block matrices, modified and stabilized 
variants and even (more recently) efficient parallel implementations 

Then m late 1970s an extensive research was under way for generalized CG 
methods that could be applied to nonsymmetiic and/or indefinite problems By 
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the eaily 1980s theie was a new plethora of new versions of geneiahzed conjugate 
giadient methods Howevei, none of these methods for geneial nonsj^mmetiic 
problems woiked in practice as well as the CG method did for symmetric positive 
definite pioblems Foi some pioblems eithei a very laige number of iteiations were 
required oi the method did not even converge In his woik Kershaw also proposed 
a generalization of the Conjugate Gradient algorithm for any non-singular coeffi- 
cient matiiv arising fiom partial differential equations Sheorey T [2001] followed 
this geneialization and developed an incomplete LU preconditioned veision of the 
conjugate gradient method that is applicable for the case of an} non-singular co- 
efficient matrix The same version of PCG has been assessed in the present study 
Some applications of this PCG algorithm can be found m Sheorey et al [2001], 
Sheoiey and Muialidhai [2001, 2003] 

Books by Axelsson [1994], Barrett et al [1993], Saad [1996], Golub et al 
[1996] and Watkins [2002] provide excellent references for iterative matrix m- 
veision in geneial and PCG in particular Pissanetzky [1984] provides detailed 
discussions of almost all popular storage formats for sparse matrix handling 

1.2 Application Problem: Czochralski Crystal 
Growth 

The growth of semiconductor single crystals is the basis for device fabrication m 
microelectronics industry Single crystals of oxides have also become important 
materials for high power laser applications These oxide crystals are found to have 
special optical, fluorescent, electrical and mechanical properties, and therefore, 
these crystals are finding new applications in a rast field of ferroelctric, pyroelec- 
tric, piezoelectric , acoustic, magnetic and opto-magnetic devices Most of these 
applications require a high level of miciostructuial and chemical perfection that 
only can be evlribited by single crystals 

Although several methods have been devised for growing single crystals, the 
Czochralski (Cz) method has virtuallv dominated the entiie pioduction of single 
crystals The popularity of this method comes from its ability to meet the strin- 
gent requirements foi puiity, doping, electrical and mechanical properties, and 
crystallographic perfection 

The basic principle of growing a single cnstal by the Czochralski method 
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IS shown m Figure 11 A crucible is initially filled with the material and then 
thermal energy is supplied by a heater surrounding the crucible This raises 
the tcmpeiature of the solid above its melting tempeiatuie A seed crystal, at 
the end of the rod, is then dipped into the melt and, after appropiiatc start- 
up procedures, slowly withdrawn from the melt To keep the feeding material 
molten, the crucible is maintained at a constant temperature above the melting 
point 

Under suitable temperature conditions re-crystalhzation in the form of a single 
crystal occurs This method bears Czochialski’s name, who in 1916 established 
it as a means to determine the ciystalhzation velocity of metals 


Pulj rate 



Figure 1 1 A schematic of the basic principle of Czochralski crystal growth pro- 
cess 

Convection in Czochralski melt is pervasive in all teriestrial growth systems 
Sources for flow include buoyancy-diivcn convection caused by the solute and 
temperature dependence of the density, surface tension gradients along fluid- 
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melt menisci, foiced convection mtioduced by the motion of solid surfaces, such as 
ciucible and crystal lotation m the Cz systems, and motion of the melt induced by 
the solidification of the mateiial These flows are impoitant causes of convection 
of heat and species and can have a dominant influence on the temperature field in 
the system and on solute mcorpoiation into the crystal A. schematic of vaiious 
tianspoit mechanisms of Czochralski process is shovn m Figuie 1 2 


Pull rate 


Crystal 

rotation 




Figuie 1 2 Schematic of vaiious transport mechanisms in a Czochralski system 


A thorough understanding of heat and mass transport is a prerequisite to the 
optimization of Czochralski ciystal giowth systems for control of crystal quality, 
as measured by the degree of ciystallogiaphic perfection of the lattice and the 
spatially uniformity of electiically active solutes m it Fluid motion in Cz-system 
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IS of concein to ciystal growers as it determines the tiansport of solute , heat, 
impuiities, etc to the giown ciystal Magnitude of stresses m the grown crystal is 
governed by the temperatuie gradients in the melt It is currently accepted that 
the tiansport properties of the melt are the piimary factors affecting the quality 
of grown crystals These properties of melt are themselves determined by a broad 
and interacting matrix of variables that include crystal lotation, crucible rotation, 
pull rate, crucible temperatuie, pressure, orientation etc Thus, state-of-the-art 
crystal growth requires algorithms that models the growth by incorporating the 
effects of all parameters and gives a set of values of these parameters as best 
possible guidelines for growing a particular crystal 

1.3 Literature Review on Crystal Growth 

Over the past tvo decades, numerical simulations of Czochralski (Cz) crystal 
growth have been performed with varying degrees of complexity Both finite 
difference and finite element methods have been used to simulate first the two 
dimensional and then the three-dimensional melt flows with many different kind of 
boundary conditions A regular cylindrical geometry has usually been considered 
m most of the cases Recently, efforts have been devoted to account for the 
ciucible bottom curvature Finite element methods can easily handle the non- 
flatness of the crucible bottom However, this task becomes difficult if a regular 
finite difference or finite volume method is used Consideration of crucible shape 
IS important because the melt in industrial Cz(Si) process can become quite 
shallow Also as the ciystal is pulled, the aspect ratio of the melt continuously 
decreases 

Other difficulties in modeling realistic Cz flovs arise from the fact that the 
positions of the melt-ciystal interface and melt meniscus are neither known a 
prion nor planar The melt-ciystal interface coincides with the freezing point 
isotherm m the sjstem (foi a pure material), and the melt meniscus results from 
a force balance between the suiface tension, pressure force and viscous effects 
To complicate this further, the flow and temperature fields which are used to 
determine the shapes and locations of interface and meniscus are often oscillatory 
m nature To overcome the difficulties associated with the time varying domains 
and irregular boundaries, different techniques have been employed by various 
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investigators 

Dei by and Blown [1986] have developed an ’’integiated hydrodynamic ther- 
mocapillaij'” (IHTCM) model based on the finite element method which allows 
the calculation of interface as well as meniscus shape The IHTCM models heat 
tiansfer bj^ conduction and diffuse, gray radiation between all components within 
the Cz puller and axisymmetnc, steady-state, laminar flow m the melt The la- 
dius of the crystal is computed simultaneously and self-consistently with the heat 
tiansfer thioughout the system and convection in the melt Kinney and Brown 
[1993] have incoipoiated k-e model for tuibulent flow into the IHTCM More re- 
cently, a global conduction-iadiation procedure based on a finite element program 
ABAQUS, and a thermal radiation model, combined with a thiee-dimensional 
NaMer-Stokes solver based on the finite element method has been developed by 
Koai et al [1994] 

In the aiea of finite differences Kopetsch [1989, 1990] extended a numerical 
procedure to include the deformable shapes of the melt-crystal interface and melt- 
ambient meniscus His numerical scheme is based on an adaptive grid generation 
that maps the irregular boundary/mterface movement and generates a coordinate 
transformation for a curvilinear domain Several papers have been published in 
this direction and many of the problems faced in the formulation have been 
resolved For example, the original scheme is suitable only for a single domain 
and steady-state problems It is however, very difficult to use this formulation for 
a multizone system since it needs an adaptive grid geneiation algorithm in the 
interface zone Also, the algorithm is not suitable for transient processes since it 
lacks the relationship between the grids at different times 

Zhang et al [1995] developed a multizone adaptive grid-generation (MAGG) 
technique foi efficient simulations of various transport phenomena on irregularly 
shaped domains and in regions with moving and/or free boundaries MAGG 
IS based on a variational optimization appioach and provides boundary-fitting 
capability and discretional y control over grid distiibution, as well a local and 
global orthogonality Based on this giid generation technique and cuivihnear 
finite volume disci etization, Prasad et al [1995] developed a high resolution 
computer model (MASTRAPP2d) to simulate crystal growth processes at low 
and high pressures This scheme allows consideration of a multiphase system for 
more than one material in one single domain which may consist of irregular and 
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moving boundaiies, interfaces and free surfaces 

Solute segiegation in the Cz systems is coupled with melt convection, species 
diffusion and the physical and chemical properties of impuiities m melt and solid 
An dccuiate simulation of segregation lemains a difficult task The major chal- 
lenge m predictrng the solute distribution in the crystal stems from the fact that 
solutal transport m melt is an intrinsically dynamic process which is strongly 
time dependent Thus dopant distribution cannot be treated as quasi-static like 
flow and teiupeiatuie fields, which do not change much m a small time interval 
since the crystal rates are generally very small 

The equilibrium segregation coefficient Lq is applicable only foi negligibly slow 
growth rates For finite or higher solidification rates, solute atoms with ko < 1 
are rejected by the advancing solid at a greater rate than they can diffuse into the 
melt Many attempts therefore have been made to obtain appropriate values of 
the effective segregation coefficient For ID (longitudinal) steady-state diffusion 
of a solute, Hurle [1993] gives 


where 


keff — 


ko 

[1 - (1 - ko) J] 


( 11 ) 




D 




(1 2 ) 


and Us IS the rotation rate for the crystal and H is a. rotation parameter 

Burton et al [1953] assumed that solute distribution varies within a thin 
layer of thickness 5 adjacent to solid-melt interface A complete mixing caused 
by convection is assumed outside this boundary They related the parameter J 
to the boundary layer thickness 5 by 


^ = -ln(l - J) (1 3) 

This approach wms criticized by Wilson [1980], wflio has demonstrated that the 
relationship should be 
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M = J (14) 

These two definitions (Equation 1 3 and Equation 1 4) approach each othei as 
J becomes very small compared with unity An evaluation of the integral m 
Equation 1 2 then yields 


5 = 1 GAc'/' (15) 
This gives the famous BPS model as 


I 2:2 (I 

[/cq + (1 T /Jo)exp(— 

where it has been assumed that Ps = pi 

In the limit when (5 — > 0, k^ff becomes equal to ko, that is, the boundary 
layer vanishes due to very strong convective flows m the melt When 5 — >■ oo, 
the effective segiegation coefficient approaches unity, that is, the mcorpoiation 
of solute m the solid is equal to that in the melt 

A detailed discussion on microsegiegation m the formalism of the boundary 
layei has been presented by Favier et al [1982] They have decoupled the convec- 
tive and diffusive material tiansport and have assumed convection to be steady 
and diffusion to be time dependent This assumption is interesting and quite 
questionable since it has been observed that unsteady convection is the cause of 
change of crystal radius 

Garandet et al [1994] demonstrated that the results of Favier’s and Wilson’s 
models arc quite close and also aie in agreement with the theoretical model of 
Hurle et al [1969] But the segregation in lateial diiection cannot be treated 
by a ID model Expeiimental results of Chung et al [1997] clearly demonstrate 
that the solute concentration is a stiong function of radial direction because of 
the strong convective effects in the melt and non-planar ciystal-melt interface 
Solute concentrations based on dynamic simulations from beginning of growth 
until the end have never been reported 
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1 4 Scope of the Present Work 

In summaiizing the pievious \\ork on Conjugate Gradient methods, it has been 
leahzed that it is an extremely effective technique for solving lineai systems when 
the coefficient matrix is symmetiic positive definite Geneialization and efficient 
implementation of Conjugate Gradient technique foi the linear systems with un- 
symmetiic coefficient matiices is still an active area of lesearch Preconditioned 
Conjugate gradient (PCG) algorithm developed by Slieoiey [2001] is one such 
geneialization and piesent lesearch is an assessment of this veisioii of the PCG 
algoiithiii The algoiithm has been examined foi the following test pioblems 

1 Two dimensional unsteady heat conduction in a region with variable ther- 
mal conductivity, 

2 Two dimensional unsteady advection-diffusion problem of flow between two 
infinite paiallel plates, 

3 Momentum, energy and solutal transpoit equations arising in the modeling 
of Czochralski crystal growth 

Effect of the preconditioning technique used i e , incomplete lower-uppei (ILU) 
preconditioning has been studied through the eigenvalue spectrum and condition 
numbers In addition, solutal tiansport has been studied with a special interest 
in understanding solute segregation in the Czochralski crystal growth process 

1.5 Thesis Organization 

The piesent thesis has been organized m the following manner 

1 Chapter 1 presents introduction, literature review and scope of the piesent 
research 

2 Chapter 2 presents mathematical background of the Conjugate Gradient 
method and then briefly summarizes various pieconditioneis 

3 Chapter 3 presents some numerical experiments with PCG 

4 Chapter 4 presents mathematical formulation for modeling transport phe- 
nomena 111 C/ochialski crystal giouth processes 
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5 Chaplei 5 gives a numerical methodology for solving the partial differential 
equations governing Gzochralski ciystal giowth 

6 Chaptei 6 caiiies conclusions of this study and scope of future woiL 



Chapter 2 

Conjugate Gradient Algorithm: 
Mathematical Background and 
Preconditioning 


In this chapter, the steepest descent method for solving linear systems is dis- 
cussed The powerful Conjugate Gradient (CG) algorithm is then presented as a 
variation on the Steepest Descent algorithm The standard conjugate gradient 
algorithm (that is applicable only when the coefficient matrix is symmetric and 
positive definite) has been generalized to the case of any nonsmgular coefficient 
matrix The general CG algorithm thus obtained is then preconditioned with 
incomplete lower-upper(LU) decomposition of the coefficient matrix Various 
preconditioneis are finally summarized 

2.1 Descent Methods 

If A IS symmetric and positive definite then the problem of solving the linear 
system Ax = b can be reformulated as a minimization problem Define a function 
J 7^" 7^ by 

J{y) = ( 21 ) 
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Then the vector that minimizes J is exactly the solution of Ax = b The can be 
seen by wilting the function J as 

Jiy) = ^V^Ay-y'^Ax 

= ^y'^Ay - y'^Ax + ^x^Ax - '^x'^Ax 

= \{v - ^{y - x ) (2 2 ) 

The rerm ^x^ Ax is independent of y, J{y) is minimized exactly when ^{y — 
x)'^A{y-x) IS minimized Since A is positive definite, this term is positive unless 
y — X = 0 Thus it takes minimum value when and only when y — x 

Descent methods solve Ax = 5 by minimizing J These are iterative methods 
Each descent method begins with an initial guess x^°^ and generates a sequence 
of iterates x^^\ such that at each step J < J 

and preferably J < J [x^'^^) In this sense we get closer to the minimum 

at each step If at some point we have Ax^^^ = 5 or nearly so, we accept 
as the solution The step from x^^^ to has two ingredients (i) choice of a 

search direction, and (ii) a line search m the chosen direction Choosing a search 
line amounts to choosing a vector that indicates the direction m which we 

will travel to get from to Once a search direction has been chosen, 

a;(^+i) will be chosen to be a point on the line | a E 7^} Thus we 

will have 

2.(fc+i) ^ ^(k) ^^p(k) (2 3) 

for some real ak The process of choosing ar^ from among all £ TZ is called 
line search Here ak is chosen is such a way that J < J One way 

to ensure this is to choose a* so that J = mmaeiiJ + ap^^'>) This 

way of choosing ak is called exact line search and it gives 


p{k)^r^k) 

pik)'^ Ap^^') 


(2 4) 


wheie = b — Ax^^^ 
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2 11 Steepest Descent 

The method of steepest descent takes and performs exact line searches 

Since = — VJ(a;^^)), the search direction is the direction of steepest descent 
of J from point 

To search m the direction of the steepest descent is a quite natural idea but 
unfortunately, it does not work particularly well Nevertheless, steepest descent is 
worth studying for at least two reasons (i) It is a good vehicle for introducing the 
idea of preconditioning (ii) Minor changes turns the steepest descent algorithm 
into the powerful conjugate-gradient method 

It IS a simple matter to program the steepest descent algorithm At each 
step approximate solution is updated using Equation 2 3 where values of ak are 
calculated using Equation 2 4 Calculation of dk requires, among all other things, 
multiplying the matrix A by the vector Cost of this operation depends on 
the degree of sparsity of A In many applications the matrix-vector product is 
the most expensive step of the algorithm and the number of these steps needs 
to be minimized Calculation of residual = b — Ax{^k.) ^-Iso seems to require 
an additional matrix-vector product Ax(^k) aiid this can be avoided by a simple 
recursion formula 


^(fc+i) ,, ^(fc) _ akAp^’^^ (2 6) 

which IS a consequence of Equation 2 3, to update the residual from one iteration 
to the next Now the matrix-vector product Ap^'^'> need not be calculated again 
as it has been already calculated as part of computation of ak A prototype of 
steepest descent algorithm is given below (Watkins [2002]) 

Algorithm 2 1 Prototype Steepest Descent Algorithm 

= 5 — Ax^^^ 

p(o) = ^(0) 

for k = 0, 1, 2, 

q{k) __ ^p{k) 

cxk = 

jj,(fc+l) _ ^{k) _|_ 
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2 12 Geometric Interpretation and Convergence 

The objective of the descent method is to minimize the function J(y) From 
Equation 2 2 one can see that J is of the form 

J{y) ^'^{y - 1 (2^) 

where x is the solution of Ax = b, and 7 is a constant Since A is symmetric 
there exist an orthogonal matrix [/ such that U'^AU is a diagonal matrix A The 
mam diagonal entries of A are eigenvalues of A, which are positive Introducing 
new coordinates z — U'^{y — x) and dropping the inessential constant 7, we see 
that minimizing J{y) is equivalent to minimizing 


J{z) — z^Kz = ^ (2 7) 

i:=l 

To get a picture of function J, consider a 2 x 2 case Now j is a function of 
two variables, so its contours or level surfaces J{zi^ ^2) = c are curves m a plane 
which have the following form 

Ai2i^ + = c, 

which are concentric ellipses centered on the origin The orthogonal coordinate 
transformation z = lF{y-x) preserves lengths and angles, so the contours of J 
are also ellipses of the same shape Such contours for the 2 x 2 case are shown in 
Figure 2 1 

The semiaxes of the ellipses are c/Ai and \/ cj A2, whose ratio is 
Eigenvalues of a positive definite matrix are the same as its singular values, so the 
spectral condition number is equal to the ratio of largest to smallest eigenvalue 
Thus the shape of the contouis depends on the condition number of A and greater 
the condition number, the more eccentric the ellipses are 

The dotted lines m Figure 2 1 represents the four steps of the steepest descent 
algorithm For a given point, the search proceeds m the direction of the steepest 
descent, which is orthogonal to the contour line The exact line search follows 
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the search line to the point at which J is minimized J decreases as long as the 
seaich line cuts through the contours The minimum occurs at the point at which 
the search line is tangent to a contoui Since the next search direction will be 
oithogonal to the contour at that point, each search direction is orthogonal to 
the previous one Thus the seaich bounces back and forth in the canyon formed 
by the function J{y) and proceeds steadily towards the minimum 

The minimum is quickly reached if A is well conditioned In the best case, 
Ai = A2, the contours will be circles, the direction of steepest descent (from any 
starting point) points directly to the center, and the exact minimum is reached 
m one iteration If A is well conditioned, contours will be nearly circular and 
the diiection of steepest descent will point close to the center, and the method 
will converge rapidly If, on the other hand, A is somewhat ill conditioned, the 
contours will be highly eccentric ellipses Piom a given point, the direction of 
steepest descent is likely to point nowheie near the center and it would not be a 
good search direction In this case the function J forms a steep, narrow canyon 
and the algorithm bounces back and foith in this canyon It takes very slow 
steps and approaches the minimum with agonizing slowness This phenomenon 
does not require extreme ill conditioning Even if the system is modestly ill 
conditioned and well worth solving from the standpoint of numerical accuracy, 
the convergence can be very slow 



2 2 The Conjugate Gradient Method 


18 


2 2 The Conjugate Gradient Method 

The steepest descent method discussed earlier is limited by its lack of memory It 
onlj' uses information about to get to and all information from earlier 

iterations is foigotten The Conjugate Gradient method is a simple variation 
on steepest descent that performs better because it has a memory Conjugate 
giadient updates the solution m each iteration by using 

+ (2 8 ) 

which IS exactl) identical to the steepest descent approach, Equation 2 3 This 
gives the following residual update 


r(*'+i) = (2 9) 

Computation of a is organized a bit differently from steepest descent, but 
the difference is cosmetic The choice of search directions, although inspired by 
the steepest descent, makes the whole difference In conjugate gradient search 
direction is updated as 


pQ H-i) _ +1) (2 10) 

where is chosen such that is A— orthogonal (conjugate) to all the previous 
search directions and this condition gives the following expression for 






( 211 ) 


In addition to the orthogonality, derivation of the above expression uses self- 
adjoint property of A with respect to the inner product, Axelsson [1994] 

After j iterations of the form both conjugate gradient 

and steepest descent algorithms give 


a;0) — Q!opf°) + 
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Thus lies 111 the space Sj spanned by the j search directions 
This IS true for both the algorithms However, they pick different search di- 
rections Intel estingly, the two different sets of the search directions span the 
same space, Sj = span {b, Ab, A^b, , A^~^b}, a Kiylov subspacc Since steepest 
descent does exact line searches, minimizes the function J ovei the j lines 
3 ,(A-) -I- ap^''\ L = 0, ,j - I The union of these lines is a tiny subset of Sj 
But, because of the conjugate seaich directions, conjugate gradient manages to 
pick out the xb) that minimizes J over the entire subspace This is called 
subspace miriimization and therefore, conjugate giadient is also called a subspace 
mimmizatioii method 

Standaid conjugate gradient algorithm can now be summarized as follows 
(Barrett et al [1993]) 

Mgorithm 2 2 Standaid Conjugate Gradient Algorithm 

1 Set A: = 0 Choose an initial approximation xq Compute = h- Ax^°'> 
and set po = tq 

2 Compute 

oik = {rk>rk)l{pk,Apk) 

_ ^(fc) 

^(K+i) — _ aiAp^^^ 

p{fc + l) _ y-ik + l) ^ 

3 Set A = /c + 1 

4 Repeat step- (2) and step- (3) until convergence in x is reached 

The method became popular due to many factors 

1 It has an optimality property ovei the relevant solution space, which usually 
means convergence to an acceptable accuracy with far fewer steps than the 
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numbei lequned foi the finite teimination piopertv — i e , it has a relatively 
high late of con\eigence 

2 The rate of comeigence can be much impioved with vaiious preconditioning 
techniques 

3 The method is paiametei fiee — i e , the usei is not lequired to estimate any 
adjustable parameters 

4 The shoit lecuiience relation makes the execution time pei iteiation and 
memoiy lequiiements acceptable 

5 The round-off erioi pioperties are acceptable 

2.3 Generalized Conjugate Gradient Method 

With standaid conjugate gradient method for solving Ax ~ h, Algorithm 2 2, 
success is guaranteed only when A is symmetric and positive definite After k 
iterations of this algorithm solution can be written in the following form 

k 

x(f') =xo + Y^ - b) (2 12) 

and the algorithm chooses ai’s to minimize \\xh — a,||, where ||^|| = (z, 

The standard conjugate gradient algorithm is modified to a generalized con- 
jugate gradient algorithm that applies to any nonsingular coefficient matrix A 
The tuck used is to solve a modified system A'^Ax = A^b and since, in this modi- 
fied system coefficient matrix is syrametiic, standard conjugate gradient becomes 
applicable Applying the standard CG algorithm to this modified system is equiv- 
alent to expanding m powers of (A^A) instead of powers of A m the light hand 
side of the Equation 2 12 and this makes us able to minimize xi — xm Euclidean 
noim For the new generalized algorithm we have the following equation 

k 

= xo + a,(A^ 4)*-' A^(Aaio - h) (2 13) 

1=1 

Generalized conjugate gradient algorithm can now be summarized as follows 
(Kershaw [1978]) 
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Algol ithm 2 3 Geneialized Conjugate Giadient Algorithm 

1 Set k = 0 Choose an initial appio\imation dq Compute = b — 
and set po = o 

2 Compute 

= {'ri,rk)l{Pk,Pi) 

2;(A.+ 1) _ 2;(/") -p 
j,(A.+ l) _ ;y(A.) _ 

p(m) ^ 

3 Set /c = /c + 1 

4 Repeat step— (2) and step— (3) until convergence in x is reached 

2.4 Preconditioned Conjugate Gradient Method 

Consider a linear system of equations, 


Ax = b (2 14) 

where A is any nonsingular sparse matrix The fact that one can always do an 
incomplete lower-upper (LU) decomposition of the coefficient matrix A is the mo- 
tivation for the selection of incomplete L{7-fdctorization as preconditioner Exact 
decomposition will be subject to fill-in and moieovei, that is almost equivalent to 
solving the linear system itself So, instead of doing complete L[/-decomposition, 
L IS defined to have non-zeio entries only within the semi bandwidth of A corre- 
sponding to the lowei triangular portion of A and similarly U corresponding to 
the upper triangular portion of A Thus the factorization is clearly approximate 
The form of L[/-decomposition algorithm in which all the diagonal elements of L 
are 1 (Crout’s Decomposition) has been used This incomplete L[/-decomposition 
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algoiithm IS summarized below 

Algorithm 2 4 Clout’s LlZ-decomposition algouthm 

1 Choose element row wise, within the band of A 

2 If Oij lies under the main diagonal, i e , i > j 

j-i 

li j j ^ k'^t ,_7 

fc=l 

3 If Cij lies on the main diagonal, i e , i = j 

1—1 

^1,1 I ^ ^ 

fc=l 

Uii — 1 

4 If a^J lies above the mam diagonal, i e , i < j 

Thus one can obtain an approximate inverse {LU)~^ for A and then, can 
rewrite the original system, Equation 2 14, in the following form 


L-KAU-^Ux) = L~^b (2 15) 

Now, the generalized conjugate giadient method, Algorithm 2 3, when applied to 
the above preconditioned system we get the following Preconditioned Conjugate 
Gradient Algorithm (Kershaw [1978]) 

Algouthm 2 5 Preconditioned Conjugate Gradient Algorithm 

1 Set A: = 0 Choose an initial appioximation xq Compute = b - 
and set po = {U'^U)~^ A'^ {LL'^)~^ r o 
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2 Compute 

= {n,{LL'^)-W,)/{p,,U'^Up,) 

3-(*. + l) _ ^(A.) ^ Q;^p{^) 

^(A + 1) ^ ^(A) _ ^^^p(A) 

= (r(^+i). (LL^)-V(^+i))/(r(^),(LL^)-ir(')) 

p(A+l) ^ ^ ^^p(A) 

3 Set A, =: A: + 1 

4 Repeat step— (2) and step— (3) until convergence in x is reached 

This preconditioned conjugate gradient is general and works foi the case of 
any nonsingular coefficient matri>^ Nevertheless, in the cases where the coeffi- 
cient matrix is quite ill-conditioned there are chances of becoming zero during 
incomplete LC/-decomposition and if so happen, algorithm will break down imme- 
diately In such a case, one can simply set to a nonzero value and go on with 
the algorithm Although this will introduce some error but Xf7-decomposition is 
approximate anyway 

2.5 Preconditioners 

As discussed earlier, the convergence rate of iterative methods depends on spectial 
propeities of the coefficient matrix Hence one may attempt to transform the 
linear system into one that is equivalent in the sense that it has the same solution, 
but has more favorable spectral properties A preconditioner is a matrix that 
affects such a transformation (Barrett et al [1993]) 

For instance, if M approximates the coefficient matiix A in some way, the 
transformed system 

il/-i 4a = M-^h (2 16) 

has the same solution as the oiiginal system Ax — 6, but the spectral properties 
of its coefficient matrix M~'^A may be more favorable 
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In devising a preconditioner, one is forced with a choice between finding a 
matrix M that appioximates 4, and for which solving a linear system a easier 
than solving one with A, or finding a matrix M that appioximates A~'^, so that 
only multiplication by M is needed The majority of available preconditioneis 
fall m the first category 

2 5 1 Cost trade-off 

Since using a pieconditioner in an iterative method incurs some extia cost, both 
initially foi the setup, and per iteration lor applying it, there is a tiade-off be- 
tween the cost of constructing and applying the pieconditioner, and the gain 
m convergence rate Certain preconditioners little or no construction phase at 
all (foi instance the SSOR preconditioner), but for others, such as incomplete 
factorizations, there can be substantial work involved 

On parallel machines there is a further trade-off between the efficacy of a 
preconditioner m the classical sense, and its parallel efficiency Most of the tra- 
ditional preconditioneis have a large sequential component 

2 5 2 Left and Right Preconditioning 

The transformation of linear system A — >■ M_iA is often not what is used in 
practice For instance, the matiix M-iA is not symmetric, so, even if A and 
M aie, the standard conjugate gradient method is immediately not applicable to 
this system In this case, a way of deriving the preconditioned conjugate gradient 
method would be to split the preconditioner as M = Mi M 2 and to transform the 
system as 

Mi-^AM2-\M2x) = Mi~^b (2 17) 

If M is symmetric and Mi = M 2 ^, it is obvious that now one has a method with 
a symmetric iteration matrix and hence, standaid conjugate algoiithm can be 
applied 

There is one more approach to preconditioning, which is easier to derive 
Again consider the system in Equation 2 17 Here Mi and M 2 aie called the left- 
and right preconditioners, respectively and one can simply apply an unpiecondi- 
tioned iterative method to this system Only two additional actions ro = Mi~^ro 
before the iterative process and = M 2 ~^Xn aftei aie necessary 
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2 5 3 Jacobi Preconditioning 

The simplest pieconditioner consists of just the diagonal of the matrix 


= 




\{i= j 


0 otherwise 


This is knoM n as point incomplete preconditioner It is possible to use this pre- 
conditionei without using any extra storage beyond that of the matrix itself 
However, division operations are usually quite costly, so in practice storage is 
allocated foi the reciprocals of the matrix diagonal 

Block Jacobi Methods 

Block versions of the Jacobi preconditioner are derived by partitioning of the 
variables If the index set S = {1, ,n} is partitioned as ^ = [J^Si with the sets 

Si mutually disjoint, then 


m- 




Oij if 2 and j are m the same index subset 
0 otherwise 


The preconditioner is now a block-diagonal matrix 

Often, natural choices for the partitioning suggest themselves 

1 In problems with multiple physical variables per node, blocks can be formed 
by grouping the equations per node 

2 In structured matrices, such as those from partial differential equations on 
regular grids, a partitioning can be caused on the physical domain Exam- 
ples are a partitioning along lines in a 2D case, or planes m 3D case 

3 On parallel computers it is natural to let the partitioning coincide with the 
division of variables over the processors 

Jacobi preconditioneis need very little storage, even m the block case, and they 
are easy to implement Additionally, on parallel computers they do not present 
any particular problems On the other hand, more sophisticated preconditioners 
usually yield a larger improvement 
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2 5 4 SSOR Preconditioning 

The SSOR (Symmetric Successive Over-Relaxation) like the Jacobi precondi- 
tioner, can be derived from the coefficient matrix without any work If the original 
symmetric matrix is decomposed as 

A — D L 

in its diagonal, lower, and upper triangular part, the SSOR matrix is defined as 

M ^ {D + L)D-\D + Vf , 

01 , parameterized by uj 

The optimal value of the w parameter, like SOR method, will reduce the number 
of iterations to a lower order Specifically, for second order elliptic problems a 
spectral condition number K,{M^opt~^A) = 0{^k{A)) is attainable, Axelsson et 
al [1984] In practice, however, the spectral information is needed to calculate 
optimal w 

The SSOR matrix is given in factored form, so this preconditioner shares many 
properties of other factorization-based methods For instance, its suitability for 
vector processors or parallel architectures depends strongly on the ordering of 
variables On the other hand, since this factorization is given a pnon, there is 
no possibility of breakdown as in incomplete factorization methods 

2 5 5 Incomplete Factorization Preconditioners 

A broad class of preconditioners is based on incomplete factoiizations of the co- 
efficient matrix The factorization is called incomplete if during the factorization 
process certain fill-tn elements, nonzero elements in the factorization in positions 
where the original matrix had a zero, have been ignored Such a preconditioner 
IS then gi\ en a factored form M = LU with L lowei and U upper triangular The 
efficacy of the pieconditionei depends on how well Tl/"’- approximates 

Unlike the Jacobi and SSOR preconditioneis, incomplete factorization precon- 
ditioners need a non-tnvial cieation stage Incomplete factorizations may break 
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down (attempted division by zero pivot) or result in indefinite matrices (negative 
pivots) even if the full factoiization of the same matrix is guaranteed to exist and 
yield a positive definite matrix 

An incomplete factoiization is guaianteed to exist for many factorization 
strategies if the original matrix is an M— matrix (z e , Mij < 0 for all z j) This 
was oiigmally proved by Meijeiink and Van der Vost [1981] In cases where the 
pivots aie zero oi negative, strategies have been proposed such as substituting 
an aibitiaiy positive number (Kershaw [1978]), or restarting the factorization on 
A + al for some positive value of a (ManteufFel [1980]) 

An impoitant consider ation for incomplete factorization preconditioners is the 
cost of the factorization process Even if the incomplete factorization exists, the 
number of operations involved in creating it is at least as much as for solving 
a system with such a coefficient matrix, so the cost may equal that of one or 
more iterations of the iterative method On parallel computers this problem is 
aggravated by the general poor parallel efficiency of the factorization 

Point Incomplete Factorizations 

The most common type of incomplete factorization is based on taking on a set S 
of matrix positions, and keeping all positions outside this set equal to zero during 
factorization The resulting factorization is incomplete in the sense that fill is 
suppressed 

The set S is usually chosen to encompass all positions (t,j) for which ^ 0 
A position that is zero in A but not so in an exact factorization is called a fill 
position, and if it is outside 5, the fill there is said to be discarded Often S is 
chosen to coincide with the set of nonzero positions in A, discarding all fill There 
are two major strategies for accepting or discarding fill-in, one structural and one 
numerical The structural strategy is that of accepting fill-in only to a certain 
level The numerical fill strategy is that of drop tolerances fill is ignored if it is 
too small, for a suitable definition of small Although this definition makes more 
sense mathematically, it is harder to implement m practice, since the amount of 
storage needed for the factorization is not easy to predict 
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Block Incomplete Factorizations 

The staitmg point for an incomplete block factorization is a partitioning of the 
matii\ Then an incomplete factorization is peifoimed using the matriv blocks 
as basic entities The most important difference with the point methods arises 
m the inversion of the pivot blocks Wheieas inverting a scaler is easily done, 
in block case two problems arise First, inverting the pivot block is likely to 
be a costly operation Second, initially the diagonal blocks of the matrix are 
likely to be spaise and one would like to maintain this stiuctuie throughout the 
factoiization Hence the need for approximation of inverses arises 

In block factorization a pivot block is generally forced to be sparse, typically 
of banded form, and that we need an approximation to its inverse that has similar 
structure Further this approximation should be easily computable, so the option 
of calculating the full inverse is generally ruled out by taking the banded part of 
it Foi example, the simplest approximation to A~^ is the diagonal matrix D of 
the reciprocals of the diagonal of A 

Banded approximations to the inverse of banded matrices have a theoretical 
justification In the context of partial differential equations the diagonal blocks 
of the coefficient matrix are usually diagonally dominant For such matrices, 
the elements of the inverse have a size that is exponentially decreasing m their 
distance from the mam diagonal 

2 5 6 Polynomial Preconditioners 

Polynomial preconditioners can be considered as the second class of precondition- 
ers where the preconditioning matrix is a direct approximation of the inverse of 
the coefficient matrix Suppose that the coefficient matrix A of the lineai system 
IS normalized to the form A = I — B, and that the spectral radius of B is less 
than one Then using the Neumann senes, one can write the inverse of A as 
A~^ = thus, an appioximation may be derived bj'' truncating this 

infinite series 

Dubois et al [1979] investigated the relationship between a basic method 
using a splitting A = M - N, and a polvnomially preconditioned method with 
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The basic lesult id that for classical methods, k steps of the polynomially pre- 
conditioned method aie exactly equivalent to kp steps of the oiiginal method 
Foi acceleiated methods, especially the Chebyshev method, the preconditioned 
iteration can impiove the numbei of iteiations by at most a factor of p 

Although there is no gain m the number of times the coefficient matiix is 
applied, polynomial preconditioning does eliminate a large fiaction of the inner 
products and update operations and there may be an overall increase m efficiency 



Chapter 3 


Numerical Experiments with 
PCG 


In this chaptei, numencal experiments with the Preconditioned Conjugate Gra- 
dient algoiithm have been presented A nonlinear unsteady heat conduction 
problem is hist considered and a study of variation in CPU time with the size 
of the problem is done for three linear system solvers, namely Gaussian Elimi- 
nation, Gauss-Seidel and Preconditioned Conjugate Gradient In the second test 
case, the suitability of Preconditioned Conjugate Gradient for solving advection- 
diffusion equation is assessed To study the effect of preconditioning, eigenvalues 
and condition numbers of the original and preconditioned matrices have been 
studied for the advection-diffusion problem 

3.1 Heat Conduction Problem 

Heie the following dimensionless equation governing unsteady heat conduction in 
a two dimensional variable conductivity region (Figure 3 1) is solved 



Thermal conductivity in this region vanes in the following mannei 


ki9) = l-P9 


(3 2) 


and the values of /? that has been considered are /? = 0 and 0 4 


3 1 Heat Conduction Problem 31 



Figure 3 1 Physical and Computational domains for the nonlmeai conduction 
pioblem 

Physical and computational domains along with initial and boundary condi- 
tions are shown in Figure 3 1 Governing Equation 3 1 is discretized using pure 
implicit scheme for the unsteady term and second order central difference for all 
other terms For a point {i,j) m the computational domain this discretization 
gives following relation for tempeiatures between two consecutive time steps 



+ (i + lesi'y'*) = 9“ 


Where s = ^ and h{= A.X = AE) is the disci etization in both X and Y 
directions 

3 11 Solution Procedure 

The solution algorithm has been used to obtain the time maichmg temperature 
field inside the computational domain 

1 Set time step p ~ 0 and initialize the tempeiature field at the time step, 
with a given initial tempeiatuie distiibution 
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2 Use the following steps to calculate temperature at (p + time step, 
by solving the nonhneai system of equations repiesented by the dis- 
cretized governing Equation 3 4 

(a) Initialize = g{p) 

(b) Evaluate k values in Equation 3 4 using = k and 

thus geneiate a linear system of equations in 

(c) Solve the system of linear equations for using a linear system 

solver 

(d) Repeat steps (b) and (c) until convergence m k and is achieved 


3 Set and p = p + 1 

4 Repeat step-2 and step-3 until steady state is reached 


3 12 Results and Discussion 

The results repoited for the nonhneai conduction problem includes 

1 Comparison of CPU time for solving the problem till steady-state for three 
solvers namely Gaussian elimination, Gauss-Seidel iterations and precondi- 
tioned conjugate gradient method 

2 Temperature contours obtained using the above linear system solvers 


Results have been presented for two variations in thermal conductivity, % e , 
/5 = 0 and 0 4 During nonhneaiitv iterations {1) the convergence criteria used is 


0l+l _ qW 


X 100 < 0 01 


foi all gild points (i) in the computational domain 

CPU time variation shown m Figure 3 2 clearlv shows that PCG is much faster 
then Gaussian elimination and Gauss Seidel Foi example, for /3 = 0 4 variation 
in thermal conductivity, problem solving with PCG is 10 — 20 times cheaper (in 
terms of CPU time) than Gauss-Seidel and almost 1900 times cheaper than Gaus- 
sian elimination Moreover, it has been generally observed that PCG becomes a 
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moie atti active choice as the number of grid points in the computational domain 
increases 

Figure 3 3 shows temperature contours inside the 2D domain foi a linear 
(/3 = 0) and nonlinear (/? =0 4) case In both the cases, temperature fields 
obtained using three different solvers are identical This verifies the uniqueness 
of inverses foi the matrices arising in this problem with respect to different matrix 
inverters used 
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Figure 3 2 Variation m CPU time with numbei of giid points in X and Y 
directions Theimal conductivity varies with tempeiature as k = l - fiO 



Figuie 3 3 Temperature contours for linear (/3 = 0) and non linear {^ = 0 4) 
conduction using three different matrix inveiteis Theimal conductivity varies 
with temperature as k = 1 ~ PO 








3 2 Advection-DifFusion Problem 


36 


3.2 Advection-Diffusion Problem 

Heie the numerical solution of the energy equation (advection-diffusion form) for 
laminar flow between two infinite paiallel plates has been consideied as the test 
problem Flow is assumed to be hydrodynamically fully developed Temperature 
of both the plates is kept constant at T = 0 for 2 : < 0, and at To = 1 for a; > 0 
Incoming fluid also has a temperature T Physical domain for the problem is 
shown m the Figure 3 4 



Figure 3 4 Laminar, hydro-dynamically fully developed flow between two infinite 
parallel plates with constant wall temperature conditions 

The velocity profile for this flow is parabolic and is given by the following 
equation 



Now with the assumptions 


1 constant thermal conductivitj of fluid, 

2 negligible viscous dissipation, 

3 negligible corapiessibility efiects 


(3 4) 


the unsteady eiieigy equation foi this flow is as follows (Bejan [1984]) 
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dT dT fd^T d‘^T\ 
dt ^ 

3 2 1 Nondimensionalization 


(3 5) 


Intioducing the following dimensionless quantities 


0 = 


T-Z 

To-Td 




t 

~ DI2 ’ 
U 


and 


where 

Tj incoming fluid temperature 
To wall temperature 
D distance between two plates 
U average a\ial velocity 

The governing Equations 3 4 and 3 5 takes the following form 


u 

U 


u 


[1 - y^] 


(3 6) 


ae dju-e) i ( bh oh \ 

dr ax Pe \ax'‘ dY^ ) ^ ' 

where 

Pe=^^ 

O: 

3 2 2 Initial and Boundary Conditions 

The physical domain, Figure 3 4 has geometiic, hydrodynamic and thermal sym- 
metry about the mid plane (the plane that contains 'll— axis and is parallel to 
plates) Therefore, the poition above this plane is considered as computational 
domain 

Now for this computational domain the non-dimensional initial and boundary 
conditions are as follows 
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Initial condition 

at r = 0, foi all A and Y 

0 = 0 

Boundary conditions 
for r > 0, 

at the inlet of the computational domain (A = 0 and 0 <Y < 1) 


0 = 0 


at the outlet of the computational domain (A'" 1 and 0 < y < 1) 



at the mid plane (A'^ > 0 and K = 0) 


de 

dY 


0 


at the upper wall (A > 0 and y = 1) 


0 = 1 


The computational domain along with dimensionless boundary conditions is shown 
below 


0 = 1 



Figure 3 5 Computational domain for the advection-diffusion problem 



3 2 Advection-DifFusion Problem 


39 


3 2 3 Disci etization 

Go\erning nondimensional equation, Equation 3 7, is disci eti/ed in the following 
manner 

Unsteady teim Pure implicit 

Advection teim QUICK (In view of minimizing eiroi due to false diffusion ) 
Conduction terms Second older cential diffeience 


NW N 


wv/< 




\v 




-® E 


® 9 

sw s 

Figure 3 6 Computational molecule foi discretization 


For this discretization, the present time-step is p and the next time-step is taken 
as p 1 A computational molecule for discretization is shown below in Figure 
3 6 The unsteady term at time p + 1, for the grid point P, can be written as 


\ / p Ar 

The advection term at point P and time p -I- 1 can be written as follows 


d{u*B) 

dX 


1 p+i 


{u*dy+^ - {u*By+^ 


J F 


AX 


(3 9) 


Using QUICK, at any time value of advection flux from the face e of the control 
volume surrounding P can be wiitten as 


where 


(u-fl). = («•«)„„, - yuRVN,(AXf + ^CC/RVUAVf 
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foi u* > 0 


and 




CURVK 


{u*9)^-2{u*e)p + {u* 0 ) 

(AA)2 


w 


CIIRVT = ^ + {'^*^)s 

" (Ay)2 


Similaily it can be written for the face 


w 


(«•«)„ = - lcURVN.(AXf + ^CURVT^{AY)'‘ 


where 


and 




CURVN,, 


{u*6)p - 2 {u*9)^ + {u*9) 


ww 


(AA)' 


GURVT^ = 




(Ay) 2 

Substituting the above values of iu*9)^ and {u*9)^ in Equation 3 9 gives 


(AA) 


. p ^ 


dx 


- (u*0)?+V] 


(3 10) 
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Disci etization for conduction terms can be wntten as 

p+i 


( d'^e 
\dx^ 


(AX)2 ^ (Ar)2 


(3 11) 


Substituting the disci etized values of various terms fiom Equations 3 8, 3 10 and 
3 11 in the governing equation 3 7, the following algebraic equation is obtained 


(SupPcc) 9^^ — (lOupPec + 24) 9^^ + [r + TupPec 4- 48 (l + 77^)] 9^^ 

- (24 - 9u*pPea) - (24i?2 _ 0P+^ - (247?^ - u* PeJ ^ 

— (ujyPec) 0^^ — (u^Pec) 6^'^y = r0p (3 12) 

Here Pec is cell Peclet number, defined as Pec = Pe(AX), R~ ^ and 

24(AAnPec 

r = — 

Ar 

3 2 4 Nusselt Number Calculation 


To addiess the basic question of heat transfer for this flow, let us consider a 
relationship between the heat flux and wall-fluid temperature diff'eience Since 
fluid temperature varies over the cross-section, AT — Tq — Tra is conventionally 
selected as the the wall-fluid temperature difference We are then interested in 
the heat transfer coefficient 


h = 


dy ) 


wall 


To -IT 


(3 13) 


Using the above heat transfer coefficient, Nusselt number can be defined as 


Nu = ^ (3 14) 

fC 

where Dh is the hydraulic diameter (D/, = 2D in this case) The value of mean 
temperature, m Equation 3 13, is given by 


T 

J. 


I-D/2 dy 

DU 


(3 15) 
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Intioducmg the previously used dimensionless quantities and using Equations 
3 13, 3 14 and 3 15, expiessions for Nusselt number aie given by the following 
two equations 


4 

Nu = (3 16) 

and 

dm^\j\{l-Y^)e]dy (3 17) 

3 2 5 Solution Procedure 

Heie Equation 3 12 represents the discretized governing equation and the initial 
and boundary conditions are as given m Section 3 2 2 Staitmg with a given initial 
condition for temperature, a marching solution m time for the entire domain has 
been found It is clear from Equation 3 12, that a linear system of equations is 
to be solved in each time-step to proceed m this way Preconditioned Conjugate 
Gradient (PCG) algorithm is used to solve this linear system at each time step In 
this manner the temperature field is obtained as a function of time until steady- 
state IS reached 

The Nusselt number is evaluated using Equations 3 16 and 3 17 Simpson’s 
1/3 rule has been used to calculate the bulk mean tempeiature given by equation 
3 17 

3 2 6 Results and Discussion 

The following results have been presented for this analysis 

1 Variation of Nusselt numbei with time at selected locations on the heated 
plate (Figuie 3 7) 

2 Steady-state variation of Nusselt number along the heated plate (Figure 
3 8) 

Both the results have been presented foi six values of the Peclet number, z e , 0 1, 
1, 10, 100, 500 and 1000 Figure 3 7 shons the variation of Nusselt number with 
time at selected A-locations along the plate foi si\ different values of Peclet 
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number Here, at any A'’— location Nusselt number decreases exponentially with 
time and then becomes constant at a steady-state v^alue This stead3''-state value 
deci eases with the distance along the heated plate This is expected since along 
the plate as the fluid gets heated, late of plate heat transfer and therefoie the 
Nusselt numbei decreases 

In Figure 3 8 both the plots depict the variation in steady-state Nusselt num- 
ber with distance along the heated plates The value of Nusselt number first 
decreases exponentially with distance and then becomes constant at a laigei dis- 
tance At this distance the flow is said to be theimally developed The value 
of Nusselt number for theimally developed flows has been reported in literature 
Foi the considered problem, analytic solutions m literature reports Nu = 7 54 
It has been observed that, when the flow is thermally developed, numerically 
calculated Nusselt number closely matches this value Moreover, at higher values 
of Peclet number it exactly matches the reported value An explanation for this 
fact is that, for this case, the analytic solution of energy equation reported in the 
literature is valid only for large values of Peclet number (Pe >> 1) 

In Figure 3 8, it is to be noted that m each plot, the upper two X — Y curves 
are slightly shifted in vertically upward direction and for these two, the values 
labeled on the K— axis are slightly inconsistent The labels on H— axis are exactly 
consistent with the bottom-most curve Labels on the A— axis are consistent 
with all the curves m all the plots 

For solving the linear system of equations, three solvers namely. Precondi- 
tioned Conjugate Gradient (PCG), point-by-point Gauss-Seidel and Gaussian 
elimination have been used Pieconditioned Conjugate Gradient (PCG) has been 
applied in two different ways In the first, incomplete low^er-upper (ILU) decompo- 
sition IS done m each time-step when the solver is used In the second, incomplete 
lower-uppei (ILU) decomposition is done only once when the solver is used for 
the very first time in the computer program For this problem. Point by Point 
Gauss-Seidel emerges as the fastest solver in terms of CPU time Performance of 
Preconditioned Conjugate Gradient (PCG), when implemented with incomplete 
LU decomposition in every time step is almost ten times slower than pomt-by- 
pomt Gauss-Seidel When it is implemented only once, it is almost four times 
slower than point-by-point Gauss-Seidel Performance of Gaussian elimination is 
awfully slow {i e , almost more than hundred times slower than point-by-Point 
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Gauss-Seidel) Moreovei, since it is a linear problem the effectiveness of PCG is 
not clearly seen The CPU time comparision is expected to show in favourable 
light foi nonlineai problems 
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Figure 3 7 Vaiiation of Nusselt numbei with time at selected A— locations, for 
different values of Peclet number Here X is the distance along the heated wall 
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3 3 Assessment of Eigenvalues and Condition 
Numbers 

3 3 1 Eigenvalues 

As discussed eailier m Section 2 5, the spread of eigenvalues of the coefficient ma- 
tii\ IS a qualitative measure of the mtiinsic difficulty in solving a Imeai system 
Thus the eigenvalue spectrum provides a qualitative basis foi comparing differ- 
ent coefficient matrices Moreover, the effect of preconditioning can be seen by 
compaiing the eigenvalue spectra of the preconditioned matrix with the original 
coefficient matrix 

3 3 2 Condition Numbers 

Another important issue m solving linear system Ax = h \s the sensitivity of 
the solution foi small perturbations in h and condition number of the coefficient 
matrix A is a measure of this sensitivity 

Consider the linear system Ax = b, where A is nonsmgular and b is nonzero 
There is a unique solution x to this system which is also nonzeio Now suppose a 
small vector 6b is added to b and consider a perturbed system Ax = b-\-5b This 
system also has a unique solution x, and let 6x denote the difference between x 
and x, so that x = x-\- 6x 

Now it can be easily shown that 

< ll-4|ll|A-'ll!|||i (318) 

Above equation provides a bound for in terms of The factor 

||A||||A"^|| IS called the condition number of A and is denoted as k{A) 

a(A) = ||A||||A-i|| (3 19) 

If the coefficient matrix is symmetric and positive definite (that is, all it’s 
eigenvalues are positive) above definition of condition number becomes 

k(A) = ^ (3 20) 

'^max 

But in geneial (for any nonsmgular coefficient matrix 4) some of the eigen- 
values will be complex and for that case, the following definition of condition 
number has been used 
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/.(A) = max |A,|/ mm |Ai| (3 21) 

I I 

Thus if «;( 4) IS not too laige, then small values of ||<5&||/1|6|| imply small values 
of ||5a;||/||xl| That is, the system is not oveily sensitive to peituibations in b 
Thus if k,(A) is not too laige, one can say that A is well condiUoned If, on the 
other hand, k(4) is large, a small value of |1(I6||/||611 does not guarantee that 
||(5T||/||a:l| will be small Since, in this case, there definitely are choices of h and 
8h foi which lesulting ||(Ia:||/||x|| is much larger than the resulting ||55j|/||&|| In 
other words, the system is potentially very sensitive to perturbations in b Thus 
if hi{A) IS large, A is said to be ill conditioned 

The best possible condition number is 1 Of course, the condition number 
(definition given by Equation 3 19) depends on the definition of norm, and a 
given matrix have a large condition number with respect to one norm and small 
condition number with respect to another As such, there is no cutoff between 
the ill-conditioned and well-conditioned matrices Any such (fuzzy) boundary 
depends upon a number of factors, including the accuracy of data used, the 
precision of floating point numbers and the amount of error one is willing to 
tolerate m the computed solution 

3 3 3 Computation of Eigenvalues using MATLAB 

Here, MATLAB function eig has been used for estimating eigenvalues This 
function produces a vector e containing the eigenvalues of A when used in the 
following format 

e = eig(A) 

For real matrices, eig(A) uses the EISPACK routines BALANC, BALBAK, 
ORTHES, and HQR2 BALANC and BALBAK balance the input matiix OR- 
THES converts a real general matrix to Hessenberg by using orthogonal similarity 
transformations Then, HQR2 finds the eigenvalues of a real upper Hessenberg 
matrix by the QR algorithm Thus the wdiole process can be summarized in 
following three points 


1 Balancing of the input matrix 
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2 Reduction of the balanced matrix to uppei Hessenberg form 

3 Calculation of eigenvalues of uppei Hessenbeig matrix using QR algorithm 
now these points can be further described as follows 

Balancing 

The sensitivity of eigenvalues to loundmg eirois during the execution of the algo- 
rithm can be reduced by the procedure of balancing The errors m the eigensystem 
found by the numerical procedure are generally propoitional to the the Euclidean 
norm of the matrix, that is, to the square root of the sum of the squares of the 
elements The idea of balancing is to use similarity tiansfoimations bj diagonal 
matrices to make rows and columns of the matrix have comparable norms, thus 
reducing the overall norm of the matrix while leaving the eigenvalues unchanged 
A symmetric matrix is already balanced 

The cost of balancing is •n? flop counts, where n is the order of the input 
matrix and the time taken by the subroutine BALANC is never more than a 
few percent of the total time requiied to find the eigenvalues It is therefore 
recommended that you always balance nonsymmetric matrices It never hurts, 
and it can substantially improve the accuracy of the eigenvalues computed for a 
badly balanced matrix 

Reduction to Hessenberg form 

To understand the need for reduction to Hessenberg foim, fiist one need to learn 
something about the QR algorithm The QR algorithm is an iterative process for 
finding eigenvalues It is based on QR decomposition, which is a direct procedure 
related to Gram-Schmidt process A single iteration of the QR algorithm is called 
a QR step oi QR iteration In a QR step for a general matrix, the cost of QR 
factorization is ^ and then the cost of matiix multiplication is 2n^ Thus 
the total cost pei QR step for a general matrix is Reduction to upper 

Hessenberg form is equally expensive and lequiies but it has to be done 
once, because upper Hessenberg form is preserved by the QR algorithm Now 
since QR iterations with the upper Hessenberg matrix aie relatively inexpensive 
and each can be done costing only O(n^) flops, QR algorithm is applied only after 
converting the input matrix to Hessenberg form 
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Algol ilhm foi Householder reduction to Hessenberg form is as follows (Watkins 

[ 2002 ]) 

Algoiithm 3 1 Householder Reduction to Hessenberg foim 
foi /c = 1 to n — 2 

vk = sign(a:)||a;||2ei 
Vk-=vil\\v,^\\^ 

■^1 +l m i m ~ ml m ~ 2Wa, '^l +l m k. m) 

^1 in,/ +1 m ~ ^+1 771 2(Ai 771 

end 

Here sign is a function such that sign(a;) returns an aiiay y of the same size 
as X, where each element of y is 

1 if the corresponding element of x is greater than zero 
0 if the corresponding element of x equals zeio 
-1 if the corresponding element of x is less than zeio 
and ei — (1, 0, , 0) 


QR algorithm 

As already described, QR algorithm is an iterative piocess for finding eigenvalues 
It is based on QR decomposition, which is a diiect procedure related to Gram- 
Sclimidt process During iterations Hessenberg foim of the matrix, is preserved 
and the diagonal elements finally converge to eigenvalues QR algorithm can be 
summarized as follows (Watkins [2002]) 

Algorithm 3 2 QR algoiithm 

If A° IS the Hessenbeig reduction of A, then 
Set A = {Q°fA°Q° 
for k = l,2 
Pick a shift ft'' 

^ _ ^kj 

Ak = + ix'^I 

end 


e g , choose 

QR factorization of A*'^ — 
Recombine factors in reverse order 
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3 3 4 Test cases for the Computation of Eigenvalues and 
Condition Numbers 

For the computation of eigenvalues and condition numbers, following two cases 
of the advection-diffusion problem have been considered 

Unsteady Nonlinear Advection-Diffusion 
Governing equation for this case is given as 


^ = i. [A (l^\ ^ (i \ 

dr dX Pe [dX dX j ^ dY ^ dY ] 

u* is assumed to be a parabolic flow profile given by 


(3 22) 


„- = 5[i-y2] 
and 


(3 23) 


k = l-(O4)0 


(3 24) 


Steady Linear Advection-DifFusion 

Governing equation for the steady and linear is a special case of Equation 3 22 


d(u*9) _ 1 / 8^9 d^9 \ 

dX ~ Pe \dX^ dY^ J 


(3 25) 


Discretization Formulation 

Governing Equations 3 22 and 3 25 has been discretized in the same way as in 
Section 3 2 3 Moreover, the initial and boundary conditions used are also exactly 
the same as given in Section 3 2 2 

This gives the following discretized equation for the unsteady-nonlinear case 

VRcfk T-r . r 3^^-pT 
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(3u^Pe,) CV - (iStipPec + 6 (j.w + dp - Ic)) C 

+ r + TiipPec + ASkp 9^^ — ^6 + 4Ap — kw^ 9upPec^ 9^ 

- [qE? {In + dkp - ks) ~ M*^Pec) C' - {}s + d.p - Aa^) - wJPec) 0^+^ 

- (uJ^PeJ Cv - («sPec) = rQ\ (3 26) 

where Pec is cell Peclet number, defined as Pec = Pe(AX), R = (AA’)/(Ay) and 

- 24(AX)Pec 
Ar 

Foi the steady-linear case it becomes 

(3u^Pec) Qww ~ (19u^>Pec + 24} + [Tu^^Pec H- 48 (l -f R^)] 6p 

- (24 - QiipPec) 9e - (24^^ - uj^Pcc) 9n - (24:R^ - ^sPec) 6s 

— (u^Pec) 9nw “ (f^sPsc) ^svv = 0 

Procedure for Calculation 

For the steady-linear case the coefficient matrix is generated using the governing 
discretized Equation 3 27 In the unsteady-nonhneai case the very first matrix 
arising just after 10 time steps is considered where the time step is taken as 0 005 
In both the cases, length of the computational domain is taken 50 and the original 
coefficient matrix has the structure shown below in Figure 3 9 
Then in both the cases preconditioned matrices are obtained from the left-right 
incomplete lower-upper (ILU) preconditioning of the original coefficient matrix 
Preconditioned matrix is given by 

P = L-UU-^ (3 28) 

Where L and U are the lower and uppei tiiangulai matrices obtained from LU 
decomposition of the original coefficient matrix A LU decomposition is incom- 
plete in the sense that complete band of the oiiginal coefficient matrix is not used 
Infact, all the possible bands (semi-bandwidths shown m Figure 3 9 by arrows) 
are have been tested for this type of incomplete LU decomposition 


J 
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Figure 3 9 Original coefficient matri\ 


Finally, these matiices are supplied to MATLAB and it calculates all the 
eigenvalues as discussed in Section 3 3 3 

3 3 5 Results and Discussion 

Following results have been presented for this analysis 

1 Variation of condition number with Peclet number for three different grid 
sizes, m the steady-linear case (Figure 3 10) 

2 Variation of condition number with Peclet number for three diffeient grid 
sizes, for the unsteady-nonlinear case (Figure 3 11) 

3 Eigenvalue spectia of the original matiices m steady-linear pioblem for 
different values of Peclet numbei (Figuie 3 12) 

4 Eigenvalue spectra of the original matrices in unsteady-nonlinear problem 
for different values of Peclet number (Figure 3 13) 

5 Eigenvalue spectra along with condition numbers foi the oiiginal and pre- 
conditioned matiices, for the steady-linear case (Figuies 3 14-3 28) 

6 Eigenvalue spectra along with condition numbers for the original and pre- 
conditioned matrices, for the unsteady-nonlinear case (Figures 3 29-3 43) 

Eigenvalue spectia showm in Figuies 3 12 and 3 13 aie foi a 101 y 11 grid 
All the other eigenvalue spectia, listed above in 5 and 6, have been piesented for 
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three difFeient giid sizes (101x11, 125x11 and 151x15) and fi\e values ofPeclet 
numbei (0 1, 1, 10, 100 and 1000) In Figures 3 14 to 3 43 there are six sets of 
spectra Spectrum (a) coiiesponds to the original coefficient matiix Spectra (b)- 
(f) coiiesponds to the different number of diagonals considered (i e , nbv l-nbw5 
respectively m Figure 3 9) during incomplete decomposition for constructing the 
preconditionei Condition numbers m above lesults have been found using the 
definition given by the Equation 3 21 Condition numbers given by Equation 
3 21 differ in magnitude from the values given by the 3 19, but show the same 
trend of variation with Peclet numbei This comparison for the original matrix 
of steady-hneai problem, for 101 x 11 grid size, is illustrated in the table below 


Peclet 

Number 

/c(A) 

Equation 3 21 

k(A) 

Equation 3 19 

0 1 

167 7373 

173 4418 

1 

160 6948 

173 3069 

10 

64 4685 

171 1982 

50 

36 3001 

152 8322 

100 

33 9591 

130 2967 

500 

35 3388 

91 9220 

800 

39 3543 

95 8429 

1000 

41 6694 

100 0026 


Thus the above table confirms the fact that the qualitative variation in con- 
dition number with Peclet numbei remains the same regardless of the definition 
of the condition number used Of course, the magnitude of the condition number 
depends upon the definition used 
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Figure 3 10 Variation of condition number with Peclet number for three different 
grid sizes in steady-linear case 



Figure 3 11 Variation of condition number with Peclet number for three different 
grid sizes in unsteady-nonlinear case 


(I 
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Figure 3 12 Eigenvalue spectra of the original matrices m steady-linear problem 
foi different values of Peclet number 



Figure 3 13 Eigenvalue spectra of the oiiginal matiices in unsteady-nonlinear 
problem foi diffeient \alues of Peclet number 



Figuie 3 14 Eigenvalue spectium and condition numbeis for Pe = 0 10, using a 
101 X 11 grid in the steady-lmeai case 
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Figure 3 15 Eigenvalue spectrum and condition numbers for Pe = Ij using a 
101 X 11 grid in the steady-lmeai case 
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Figure 3 16 Eigenvalue spectrum and condition numbers foi Pe = 10, using a 
101 X 11 grid m the steady-linear case 
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CN = 1 1479 



CN= 1 0000 



Figuie 3 17 Eigenvalue spectium and condition numbers for Pe — 100, using a 
101 X 11 grid in the steady-linear case 
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Figure 3 18 Eigenvalue spectrum and condition numbeis for Pe = 1000, using a 
101 X 11 grid m the steady-linear case 
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Figure 3 19 Eigenvalue spectium and condition numbers foi Pe = 0 10, using a 
125 X 11 grid m the steady-lineai case 
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Figure 3 24 Eigenvalue spectrum and condition numbeis for Pe = 0 10, using a 
151 X 15 gild in the steady-linear case 
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Figure 3 26 Eigenvalue spectium and condition numbers foi Pe = 10, using 
151 X 15 gild in the steady-lmeai case 
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Figure 3 31 Eigenvalue spectium and condition numbeis for Pe = 10, using a 
101 X 11 grid in the unsteady-nonlinear case 














Figuie 3 32 Eigenvalue spectrum and condition numbers for Pe = 100, using a 
101 X 11 grid in the unsteady-nonlineai case 
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Figure 3 33 Eigenvalue spectium and condition numbers for Pe = 1000, using a 
101 X 11 grid in the unsteady-nonlmeai case 
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Figure 3 34 Eigenvalue spectium and condition numbers foi Pe = 0 10, using a 
125 X 11 grid in the unsteady-nonlinear case 







3 3 Assessment of Eigenvalues and Condition Numbers 


^8 



Figure 3 35 Eigenvalue spectrum and condition numbers for Pe = 1, using a 
125 X 11 grid in the unsteady-nonlmeai case 














Figure 3 36 Eigenvalue spectrum and condition numbers for Pe — 10, using a 
125 X 11 grid m the unsteady-nonlineai case 
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Figure 3 37 Eigenvalue spectium and condition numbers for Pe = 100, using a 
125 X 11 grid in the unsteady-nonlineai case 
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Figure 3 38 Eigenvalue spectrum and condition numbeis for Pe = 1000, using 
125 X 11 grid in the unsteady-nonhneai case 
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Figure 3 39 Eigenvalue spectrum and condition numbers for Pe — 0 10, using a 
151 X 15 grid in the unstead 3 f-nonlinear case 









Figure 3 40 Eigenvalue spectrum and condition numbeis for Pe = 1, using a 
151 X 15 grid in the unsteadv-nonliiiear case 
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Figure 3 41 Eigenvalue spectrum and condition numbeis for Pe = 10, using a 
151 X 15 grid in the unsteady-nonlineai case 












Figure 3 42 Eigenvalue spectium and condition numbeis for Pe = 100, using a 
151 X 15 grid in the unsteady-nonlineai case 
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Fiom the results piesented in this study, the following observations can be 
lecoided 

1 As shown in Figuie 3 10, condition-number for the original matrix in the 
steady-lineai pioblem, initially decreases with an increase in Peclet number 
(for Peclet number 0 10 to 100) Latei, it staits increasing with an in- 
crease in Peclet number (for Peclet number > 500) In between, for Peclet 
number 100 to 500, a transition from decreasing to an increasing condition- 
number trend takes place Condition-number for the original matrix in the 
unsteady-nonlinear problem (Figure 3 11), decreases almost asymptotically 
with Peclet number and then, finally becomes constant with a magnitude 
one 

Here, the mathematical nature of the partial differential equation being 
solved depends on Peclet number As the Peclet number increases the 
problem changes from elliptic to hyperbolic Since we are solving it as 
an elliptic problem for all values of Peclet number, obtaining the solution 
becomes difficult foi higher values of Peclet numbers where the problem is 
actually hyperbolic Rate of convergence in our solution procedure validates 
this point 

Condition numbers represent the mathematical complexity m solving the 
problem In Figure 3 10 and Figure 3 11, when the condition is decreasing 
with increasing Peclet number, mathematically the problem is changing 
from an elliptic to relatively simple hyperbolic problem But of course, as 
observed from the convergence experience, practically solving the problem 
IS becoming increasingly difficult with increasing Peclet number 

2 Eigenvalue spectra shown m Figure 3 12 and Figure 3 12 clearly show an 
increase in the scatter of eigenvalues with increase m Peclet number There- 
fore, the relatively poor convergence of iterative methods for highei values 
of Peclet numbers is explained 

3 In both the problems consideied, condition numbei increases wnth a de- 
crease m gud spacmg foi a given lalue of Peclet number (Figure 3 10 and 
Figure 3 11) This effect is relatively less prominent at gieater values of 
Peclet number especially toi the unsteady-nonlinear pioblem It han been 
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fuithei obseived that, although the eigenvalue spectium for these matrices 
remains appio\imately same on finei giids, a few eigenvalues move closer 
and closer to zero This effect on finer giids results in an increased condition- 
numbei and gnes use to pooi convergence of iteiatue sobers 

4 Effect of /LtZ-pieconditioning is well seen foi both the problems consid- 
ered, Figures 3 14 to 3 28 foi the steady-lmeai problem and fiom Figures 
3 29 to 3 43 foi the unsteady-nonlinear problem As the width of the band 
considered during LU decomposition is increased, the eigenvalue spectium 
improves (i e , the spread of eigenvalues decreases) and condition number 
also decreases This observation can be explained from the fact that as 
the number no diagonals considered during LU decomposition increases, 
eiror introduced during LU decomposition decreases and the decomposi- 
tion becomes more and more exact Therefore, the preconditioned matrix 
L~^AU~^ becomes more closer to the identity matrix I 



Chapter 4 


Mathematical Formulation for 
Modeling Transport Phenomena 
in Czochralski Crystal Growth 
Processes 


Czochralski crystal growth is a physically complex and mathematically rich pro- 
cess The quality of grown ciystals depends on many macroscopic and microscopic 
physical effects (Figure 4 1) Formulating a mathematical model to incorporate 
all these effects is not an easy task and judicious choices need to be made m this 
regard Mathematical model presented here allows the consideration of as many 
effects as the problem demands and is believed ciucial for the modeling process 
This model diaws heavily from the earlier works of Zhang et al [1995], Prasad 
et al [1997] and Zou et al [1997] and ongoing Ph D thesis of Banerjee [2004] 

Figure 4 2 shows a schematic of the right-half of the considered axisymmetric 
Cz-system The domain consist of various mateiials in different phases with 
significantly diffeient thermophysical and transport properties (Table 4 4) The 
Cz-model presented here focuses on the growth of Nd YAG (Neodymium doped 
Yttrium Aluminum Garnet) single crystals Nd\AG falls in the category of 
oxide crystals, and since, the melt phase for the giow th of these crystals does not 
contain any volatile component, need foi an encapsulant layer ovei the melt gets 

eliminated 

An introduction to the basic growth process and various tiaiisport mechanisms 
that take place inside a Cz-system is given in Section 1 2 of Chapter 1 A list of 
vaiious important process parameters is given in Table 4 3 
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Figme 4 1 Transpoit mechanism and mhomogeneities associated with crystal 
giowth processes 

4.1 Equations governing flow, heat transfer and 
solute transport 

The assumptions that go into developing the mathematical model for flow, heat 
transfei and solute transport inside the considered axisymmetric s;>stem are as 
follows 

1 The fluid flow is two-dimensional, laminar and incompressible 

2 The fluids are Newtonian 

3 Thermophysical properties aie constant and uniform in various phases 

4 Boussmesq approximation for buoyancy-driven convection is applicable 

5 Viscous dissipation is negligible 

6 Free surface defoi matron is negligible, i e , free surface is assumed to be 
flat 
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7 Change m surface tension of the liquid phase ith temperature is negligible, 
i e , Maragoni convection is absent 

8 Liquid phase is treated as a dilute solution of the solute in the melt 

9 Segregation coefficient is independent of the local growth rate of the ciystal 
and other operating conditions 

10 No diffusion takes place m the solid phase 



Figure 4 2 Schematic of the right-half of the axisymmetnc Cz-system considered, 
showing various zones, interfaces and the free surface 

The scales that are used to nondimensionahze the governing equations are 
given in Table 4 1 Where ip refers to the dimensionless thermophysical proper- 
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ties, such as, p, p. P, /3 C, and k Subscript o refer to the reference phase, t e , 
the melt phase at freezing temperature Ty and supeiscnpt represents dimen- 
sional quantities (See nomenclatuie) Foi buovanc\ -induced flows, Uo = 
can be used for non-dimensionalization leading to Re = 1 Under the foregoing 
assumptions, the equations for mass, momentum, energy and solute transport, as 
used by Zhang ei al [1996], aie as follows 

Mass 


^ -f V {p^u^) = 0 (4 1) 

Momentum 

+ v teu,u,) = j(^V (ftVu,)- VP + A^e (4 2) 

+ v (AU.u,) = iv (AVu.)-VP (4 3) 

at Re 

Energy 

+ v (AC,.u.e) = j^v (fc.v9) (4 4) 

Solute transport 

^ + {u VC) = ^V2C (4 5) 

ut be 


In the above equations, u is the velocity vector with u and v its components 
in the X (or z) and y (or r) directions, respectively 
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Table 4 1 Scales foi Nondimensionahzation 


Variable 

Nondimensionahzation 

length 

ii,y) = {x\i/)/L 

time 

t = t*l{L/Ua) 

velocity 

{u,v) = {u*,v*)/bo 

tempeiature 

e = {T-Tf)/{n-Tf) 

pressure 

P= (p*- pgr* - Pa)/{pUo-) 

concentiation 

c = c*/c;,j 

theimophysical piopeities 

A = 


4.2 Generalized conservative form of the trans- 
port equations 

The governing equations piesented in Section 4 1 can be further put in the fol- 
lowing geneial mathematical form, Zhang et al [1995] 


d{r'^p(j)) 

dt 


+ (4 6) 


where 0 is the generalized variable, S is the volumetric source, and P is the 
diffusion coefficient The index q is set to zeio if Cartesian cooidmates aie used 
and IS unity for polar coordinates Equation 4 6 can be used foi the entire multi- 
phase, multicomponent domain with the piovision to account foi local pioperties 
and abrupt changes m tiansport properties across the zone boundaries and their 
possible movements The values for (p, P and S foi continuitv momentum, energv 
and concentration equations for an axisv mmetric geometry aie given m Table 4 2 
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Table 4 2 

Vaiiables of Equation 4 6 in an axisymmetric geometr^ 

Equations 

Variable 4) 

Diffusi\it\ r 

Souice term S 

continuity 

1 

0 

0 

^-momentum 

V 

Ih 

-dPid. + GiAa^ 

r-momentum 

V 

fil 

-dPjdr - piU/P + piw'^lr 

? m-momentum ? w 

Pi 

— (2/Ii/r)5(ni;)/5r 

0-energy 

Cpi9 

li/{Cp,Y>x) 

0 

concentration 

C 

1/Sc 

0 


4.3 Initial and boundary conditions 

The governing equations given in above sections require both initial and boundary 
conditions for a complete solution Detail of these conditions is summarized 
below 

4 3 1 Initial Conditions 

Initial conditions for the set of governing equations depends on the solution strat- 
egy adopted to get the complete solution In case, a steady-state solution is sought 
any initial field can be taken as the initial condition However, for transient anal- 
ysis exact initial conditions need to be used These aspects are discussed in 
Chapter 5 along with the solution strategy 

4 3 2 Boundary conditions 

Boundary conditions at various boundaries and intei faces in a Cz-system aie as 
follows 


4 3 Initial and boundary conditions 


95 


Bottom and side wall of the crucible 

Here foi the velocities, no-slip boundaiy condition is applicable Which can be 
wiitten as 


u = v = 0, w = Rcc? (4 7) 

and for the energy equation we have specified temperature at these walls and 
that gives the following nondimensional boundary condition 


9^1 (4 8) 

For the concentiation equation free gradient boundaiy condition is applicable 

(VC) n = 0 (4 9) 

Top and side wall of the enclosure 

Again foi the velocities, no-slip boundary condition is applicable That is 


u = V = w — 0 (4 10) 

and for the energy equation we ha\e specified temperatuie variation at these 
walls 


9 — specified variation 


(4 11) 


For the concentration equation fiee gradient boundaiy condition is applicable 


(VC) n = 0 


(4 12 ) 


Axis of symmetry 

For the portion of the axis that is inside melt following conditions aie applicable 
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de 

dr 


= 0 


(4 14) 



Solid poition of the axis of symmetry need only the following 


(4 15) 


^ = 0 (4 10) 

ai 

Melt-crystal interface 

Heie for the velocities, no-slip boundar-y condition is applicable Which can be 
written as 


ti = r; = 0, w = Re^r (4 17) 

and for the energy equation we have the fusion temperature T/ at these walls and 
that gives the following nondimensional boundary condition 


^ = 0 (4 18) 

At the melt-crystal interface, solute conservation gives the following boundary 
condition 


-^(VC) n = ^(1 - A:o)Vpuii n (4 19) 

oC 

wheie ko is the segregation coefficient and Vpuu is the nondimensional pull late 
Above boundary condition has been deiived b} equating the segregation flux at 
the interface with the diffusion flux just below it in the melt Detailed discussion 
of the same can be found in Section 5 8 of Chapter 5 
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Free surface 

Heie with the flat inteiface assumption, the kinematic and zero shear stress 
boundaiy conditions at the free surface gives the following 


dv 

11 = 0 , — = 0 

oz 


foi the sivirl velocity 


dw 

dz 

A flux balance at the free-suiface gives 

dT 


-k 


dz 


. dT 
^ dz 


+ ea,(T' - T4) 


which on nondimensionalization gives the following 


(4 20) 


(4 21) 


(4 22) 


-k^-^ 

^dz 


‘ dz 


+ Blr,;g(0/p — doo) 


and for the concentration equation again the following is applicable 


(4 23) 


dr 


= 0 


Top and side wall of the crystal 

Here we no slip foi the all the \elocity components 

U = Us, V =0, w = RQsI 


and for temperature, the following flux balance 

+ Bly s(0s ~ ^oo) 


dB 

-^Yz 


~ ^dz 


(4 24) 


(4 25) 


(4 26) 


Ik 
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4 4 Process parameters and material properties 

Tianspoit phenomena inside the and tlieicfoi^' the cijstal quality 

depends on the effect of a numbei of s\&tcm paiametcm A. list of the various 
important paiameteis inside the C/-system is gi\en bclo\* in Tabic 4 2 Typical 
yalucs of mateiial properties for Nd YA.G single citstals is giyen in Table 4 3 

Table 4 3 Various Piocess paiameteis m a C/-sy stern 


Pull rate 

Initial melt height 
Enclosuie diametei 
Crucible temperature 
Ambient temperature 


Crystal diameter 
Initial ci\stal height 
Enclosure height 
Enclosure temperature yariation 
Gas pressure 


Crystal rotation 


Crucible rotation 
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Table 4 4 Theimophvsical piopeities of Nd ^ AlG single cr^ stals m here subscripts 
6, I, q and c denote solid-ciystal, melt gas and ciucible lespectneh 


Description 


Symbol Value 


Density (solid) 

Density (liquid) 

Density (gas) 

Thermal conductivity (solid) 
Theimal conductivity (liquid) 
Thermal conductivity (gas) 
Melting point 
Melt expansivity 
Viscosity (liquid) 

Viscosity (gas) 

Heat Capacity (solid) 

Heat Capacity (liquid) 

Heat Capacity (gas) 

Heat of fusion 
Emissivity (solid) 

Emissivity (liquid) 
Emissivity (gas) 


p, 4300 Kg/in^ 

Pi 3600 Kg/m^ 

Pg 0 1602 Kg/nD 
8\V/(mK) 
ki 4 W/(mK) 

kg 0 139 \V/(mK) 

Tf 2243 K 

p 18 xl 0 “®K-^ 

m 4 68x10"^ Kg/ (ms) 

Pg 6 93x10'^ Kg/ (ms) 

c, 800 J/(KgK) 

Cl 800 J/(KgK) 

c, 1419 J/(KgK) 

AH/ 4 55xlOM/Kg 



Chapter 5 

Numerical solution of the partial 
differential equations governing 
Czochralski crystal growth 


In this chaptei the numerical methodology for solving the generalized conservative 
form of the governing equations, Equation 4 6, is presented (Baneijee [2004], 
Zhang et al [1995]) Then, the overall solution approach is discussed along with 
the details of the calculations for crystal pull rate and solute segregation 

5.1 Coordinate transformation 

A physical problem composed of a domain of irregular shape is difficult to solve 
in a regular orthogonal coordinate system such as Cartesian owing to the fact 
that the boundaries of such domains may not conform to the regular coordinate 
axes Hence an accurate representation of the system geometry becomes difficult 
It is m this context that the concept of boundaiy fitted (body- fitted) coordinate 
system has evolved When one use such a coordinate transformation, an irregu- 
lar geometry in a physical space (such as Cartesian) is transfoimed to a regular 
geometry in a computational space v herein boundaries of the transformed ge- 
ometry conform to the coordinate axes and then calculations aie performed m 
the computational space This transformation, which depends on the shape of 
the physical domain of interest, does not necessaiily guarantee orthogonality of 
the grid lines in the transformed plane In fact in most cases, the transformation 
results m a generalized nonorthogonal curvilmeai system in which the problem 
is to be solved 
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Introduction of nonoithogonahty is generally detrimentdl as the transformed 
governing equations contain extra terms to account foi non-orthogonality of the 
cooidmate system However, the fact that an accurate representation of the phys- 
ical domain is possible by this technique largely outweighs the abo^e mentioned 
limitations 

The conseivation equation 4 6 in a generalized cooidmate svstem (^,?]) can 
be written as 


Jf - ft JJ + IjKft - ftftl = r-'JS^ (5 1) 

wheie 


- -r 


r fdcp 




and 


Jr) P^Tj^ 

Othei coefficients are 

= r^h^hn'^/J 


r_ 

hj) 


dvj 


ar, = r'^hrjh^'^/J, 


= r"A/i,/ J, Pr) - r’^XhJJ, 


A = X^Xr, + y^Vr), J - X^Vn - y^Xr, 

It IS important to note some of the properties of coefficient of o- and /3 Both 
these coefficients have units of area, and while a s are alwa\s positive, /3’s may 
be positive, negative or zero The absolute magnitude of 6^ is always less than 
q;^, and if the coordinates (^, rp are orthogonal, then = 0 and = hr, These 
coefficients can be conceptually regarded as areas In view of their properties, q s 
are referred to as the primary areas and d’s as secondary aieas The subscripts 
^ and 77 in these coefficients denote their coirespondence to a ^ = constant and 
r] = constant lines respectively With this mteipietation the flux across a finite 
volume face may also be thought of as composed of two components 
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5.2 Grid generation 

As IS e\ident, Cz gl 0 ^^th &\stems die composed of zones \ itli \astl} diffeient 
thermophysical piopeities and inteinalh mos'ing mteifacps and fiee suifaces To 
peifoim accuiate simulations foi these puiposcs the giid s\stem in different 
zones need to be generated in such a way that the interfaces separating these 
zones are preserved and coincide with some grid lines, peimitting grid points to 
move only along these interfaces Here for grid geneiation, MA.GG (Multizone 
Adaptive Giid Generation) scheme developed bv Zhang ct al [1995] has been 
used The scheme has been de\ eloped based on the variational method and 
minimizes an integral function w Inch is a measure of grid chai actenstics namelv 
the smoothness, orthogonality, weighted cell area and inertia of the grids The 
scheme allow^s grids to mo\e adpatively as the solutions piogiess domains change, 
or both Moreover, the giids aie concentrated in the regions of laige variations 
m field variables using appropriate weighting functions 

5 3 Finite volume discretization 



Figure 5 1 Cuivilinear finite- Volume grid arrangement 
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The gilds used aie a structuied trapezoidal mesh, as shown in Figure 5 1 For 
a typical point P, the generalized consenation equation Equation 5 1 can be 
discietized in the following manner 

- (afJt - iSM.Ar, 

+(a,J„ - - (a,J, - = (.V.SjpAfAr, 

A.bove equation repiesents an o\eiall consei\ation of d m a finite volume in terms 
of its fluxes at the auxiliary nodes and is lefeiied to as the flux discretization 
equation The souice teim S is the original source teim plus the terms consisting 
of the gild velocity components, 


= + (53) 

07] 

where the grid velocities are defined as 

_dydx (5 4) 

~ drj dt dt dt 

Although in the calculation grids velocity components have been neglected 

and 5 = 5 IS used 

Now to obtain the finite difference form of Equation 5 2 in terms of (p at the 
centers of the finite volumes. Power Law' scheme of Patankar [1980] is used The 
fluxes using this scheme can be represented as follows 


— Fgpp + [De-^{\P^e\) '^ [\ Fe,0|]](^p Pe) 

^ ^^o(t>P + [DwM\Peu>\) + [| 0 |]](dw - 

= FnpP + [DnA{\Pen\) + [I -Fn 0 OKdn - M 
{ar^Jr,)s = F.pp + [Ds 4(|Pe.l) + [I 0, F, \]]{(ps - <Pp) 

where 




Fe = {Pn)e^ = 


Zk 

D. 


Lt / c— 


A/: 


(3^/rp +3^/rc) 


(5 6) 
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and the function A{\Pee\) is defined bv Pataiil ai [lOSO] 

A(|Pee|) = max [O 1 - 0 l|Fcl^] (5 7) 

The othei coefhcients along faces w, n, s are similar to Equation 5 G Substi- 
tuting these fluxes fiom Equation 5 5 into Equation 5 2 \^e obtain the following 
discretized governing equation in the rj coordinates 

(Xp(j)p = O-E^E T T O-N^N T OjS<t>S ^ (5 8) 

where 

aE = peA(lPeel) + [I -Fe, 0 |]] (0P - 0ir), (5 9) 

avK = [D^Ai\Pe^\) + [1 0, |]] - 4>p), (5 10) 

aj, = [DnA{\Per^\) + [1 -F„,0 1]] {(I>p - (5 H) 

as = [DsA{\Pes\) + [| 0, F, 1]] (ds - ^p), (5 12) 

„o _ (5 13) 

At ^ 

ap = aE P aw + a^^ + as + a°p — SpJA^Ar] (5 14) 

The remaining terms from the substitution aie included in b In order to identify 
the origin of these terms 5, is expressed as 

where 

(5 16) 


bs = ofp(f>p + ^cJ 
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and 


^no — ~ f 

-fFe - F, jppAr/ ~(Fn-F )opAe 

In the above discietued formulation the source tcim is lineaii^ed as = Sc + 
Sp(l>p Terms corresponding to the grid velocit} components (namely, the late of 
change of grid position v ith time) hd\e been neglected 

5 4 Treatment of pressure-velocity coupling 

Staggered girds have commonly been used for transport problems and are con- 
sidered to be more appropiiate foi pressure-^eloclt\ coupling Howe\er there are 
many advantages of using a non-staggered grid svstem for the present types of 
problem For example, the imposition of intei facial boundary conditions m a non- 
staggered grid IS more appropriate The only diavback of this approach is that 
it requires a higher order interpolation to calculate fluxes and handle pressure- 
velocity coupling The solution algorithm used for fluid flow calculations in the 
general cunilmeai cooidinate system is basically similar to the SIMPLER al- 
gorithm Patankar [1980], which consist of solving a pressure pressure equation 
to obtain the pressure fleld and solving a pressure-correction equation to correct 
the predicted velocities Howexer the scheme is more complicated because the 
velocity directions change continualh along the cooidinate lines 

5 5 Overall solution scheme for solving the com- 
plete set of equations 

The overall solution scheme foi solving the complete set of governing equations 
for Czochralski crystal growth is as follows The growth process is simulated b\ 
a senes of quasi-static solutions for flow and temperature fields and the solute 
transport equation is solved m a time niai clung wax for the whole growdh period 
This approach is based on the obsenatioii that the variation of flow and tem- 
perature fields IS fairly slow' with the growth process and theiefoie, a quasi-static 
solution shall catch the basic features of the process On the other hand, the 
solute partitioning at the melt-crystal interface and its transport m the melt aic 
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intrinsically dyndiiiic processes which are strongly time dependent Theiefore, 
a transient solution need to be obtained to tiacf the solute partitioning and its 
ledistribution in the melt duiing giovth The quasi-static solution is obtained 
foi the flow fleld at a selected ciystal height and then, based on this flov held the 
concentiation equation is solved in a time marching waj pei forming a senes 
of such quasi-static calculations, a period of growth piocess can be simulated and 
influences of vaiious related factors can be in\estigdted 

5.6 Numerical treatment of moving interfaces 



Figure 5 2 Movement of melt-crystal interface and free-surface 


Along with the solution of the fleld equations (mass, momentum and energy), 
motion of the melt-ciystal interface and the free-surface need to be determined 
Solidification of the melt is modeled as a pure substance with a fixed fusion 
temperature T/ This implies that the solid and liquid phases are separated by a 
sharp interface. .{z,r,t) = = where H is the dimensionless height 

of the melt-crystal interface, Figure 5 2 Energy balance at the inteiface (Stefan 
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condition) defines its position and motion and can be written as 


dt 


-U, 


putl 


_ Ste; 7 - 89 s d9i \ 

“ Pp 


(5 18) 


where Upuu is dimensionless crj^stal pull \elocitj in z-direction 

Pull velocity calculation is based on the assumption that the cn stal and 
melt aie not separated at the tri-junction Uiow is the dimensionless free-suiface 
velocity m —e^ direction and Usoi is the dimensionless solidification lelocit'v in 
— n diiection, Figure 5 2 Now, melt-cr\stal attachment assumption gi\es the 
following at the tii-j unction 


[Usoi{e. ii)]|^ = + {Up^a)\^ (5 19) 


Overall mass conseivation for solidification gives the following condition 


Uio^\^^Rx'^p,[Usoi{e. n)]|^ (5 20) 

where Rr is the radius ration (crystal-to-crucible) Solving above three equations 
gives the following expiession for the crystal pull velocity 


Upiiii — (1 psRr ) 


/ 




__i _ ^'] (e 
dn dn ) ^ 


n) 


(5 21) 


5 7 Solute segregation 

During the ciystal growth process, at the melt-cijstal mteiface, there is a ten- 
dency foi some of the solutes to lemaiii in the melt and foi the othcis to prefer 
the solid This phenomenon causes the concentration of the solutes to accumu- 
late or decrease and is teimed as solute segregation Howe\ei 5 convection iii the 
melt produces mixing and alters the chaiactenstics of the diflFusion layer adjacent 
to* the mteiface The spatial stiucture and intensity of the flow determine the 
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axial and lateral (perpendicular to the growth direction) profiles of the solute 
concentration in the ciystal 

Solute concentiation of the solid phase is different fiom the liquid phase due 
to segiegation The equilibrium segregation during the solidification of a binary 
system can be deteiinined fiom the coriespondmg phase diagram In many cases, 
for low solute concentration the solidus and hquidus curves on a phase diagram 
can be consideied as straight lines in the vicinity of melting temperature This 
implies that the ratio of the solubility of a solute in the solid C'„ to that in the 
melt Cl, remains constant over a concentration range This ratio is referied as 
the equilibrium segregation coefficient 


Cl 


(5 22) 


Segregation coefficient for the considered Nd YAG system is ^ 0 2, Ln et al 

[ 2000 ] 


Crystal 


Gas 

Segregation Flux 

Figure 5 3 Solute Segregation Through the Melt-Crystal Interface 


Melt 


5.8 Solution of the solute transport equation 

During the transient solution of the solute transpoit equation, initial dimension- 
less concentration of Nd in the melt is taken as initial condition and boundary 
conditions used are as given in Section 4 3 of Chapter 4 Further elaboration on 
the boundary condition at the melt-crystal interface is as follows 

Actually here the argument is that if a melt of solute concentration C* so- 
lidifies then according to the definition of equilibnum segregation coefficient, a 
solid ciystal of concentiation k^C* will be formed Since ko is less than one, some 
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of the solute must come back in the melt through melt-crystal interface This 
phenomenon is known as solute segregation and segregation flux at melt-crystal 
inteiface can be calculated using Pick’s law of diffusion Thus in dimensional 
form we have the following equation 


-D{VC*) n=:C*(l-ko)V;„H n (5 23) 

and this on non-dimensionahzation gives the follow mg 

-1(V6') n = C(1 - ko)Vpu„ n (5 24) 

wheie ko is the segregation coefficient and Vpuii is the nondimensional pull rate 
Two major problems are encountered during the solution of the solute trans- 
port equation The first is the implementation of the above segregation boundary 
condition at the melt-ciystal interface In a solidification process, solute is re- 
jected into the melt if the segregation coefficient is less than unity This increases 
the solute concentration in the melt near the interface To address this issue, 
the appioach of introducing an internal source term to simulate the segregation 
phenomenon has been followed The solute is rejected as a flux, (1 - k)VpuiiC 
and we have converted this flux to a volumetric source in the fiist control volume 
cells in the melt near the interface, such that for any of such cells 

solute input rate due to segregation = solute generation rate within the 

cell due to the source term 


This gives source term strength within a cell as 


Scell — C[1 — ko)Vpull X 

where A is the non-dimensional area, vol is non-dimensional volume of the cell 
and Vp^a is the vertical component of pull rate, Vp,.„ The possible error that 
IS introduced because of this approach is minimal m this model since the grid 
height near the interface is veiy small due to adaptive natuie of the grid and its 

clustering 
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Another issue m sohing the mass transport equation is to control the total 
mass in the system During the growth, the melt domain continuously drops 
as moie and moie melt is solidified into the cnstal This is ticated b\ simply 
legeneiatmg a new giid sjstem in the changed domain Therefore it cannot 
guaiantee the total mass conservation after each new gnd system is generated 
In addition, concentiation interpolation is lequired to stait calculations on the 
new grid Although a very little imbalance is introduced each time, a large 
effect in mass deficit may be produced ovei a long growth period because of the 
accumulating effect To avoid this difficult} the approach of conserMiig solute 
mass by pioportionally coriectmg the solute concentiations in each grid has been 
followed 


5.9 Code Validation and Grid Independence 

The computer code has been extensively validated, for example against the nu- 
merical results of Kobayashi [1978] and Prasad et al [1997] These code vali- 
dation and grid independence studies have been reported m the progress report 
of BRNS project entitled, “Mathematical Modeling of the CVD and Czochralski 
Crystal Growth Systems, Role of Magnetic Fields, and their Effects on Thermo- 
mechanical Properties” Giid independent solutions were obtained on a 180x60 
mesh 


5.10 Results and Discussion 

Using the solution scheme discussed in Sections 5 5 and 5 8 following numerical 
results have been obtained for the growth of Nd YAG single crystals 

1 Quasi-static temperature and flow fields in the melt at three melt heights 
and crystal rotations (Figure 5 7-Figure 5 9) 

2 Variation in pull-velocity of a continuously growing crystal for different melt 
heights and crystal rotations (Figuie 5 10) 

3 Variation in concentiation field ... the melt with lime for different ciystal 
rotations (Figure 5 11-Figuie 5 16) 
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4 Radial vauation m solute concentiation at different times foi different crys- 
tal rotations (Figure 5 17-Figure 5 18) 

Moieoiei, the following va'ues of vaiious important dimensionless paiameters 
have been used 
Crucible ladius = 1, 

Enclosure radius = 1, 

Crystal radius = 0 5, 

Initial crystal-height = 0 1, 

Total height ot the domain in axial direction = 3 0, 

Initial melt-height = 1, 

Crucible rotation = 0, 

Crucible temperature = 1, 

Ambient temperature = -5 06 

The enclosure wall temperature drops linearly fiom crucible temperature to am- 
bient temperature towards top 

The results have been obtained using Preconditioned Conjugate Gradient 
(PCG) as the linear system solver Results obtained were found to be exactly 
identical to the Gauss-Seidel as well as the Line-by-hne TDMA iterative solvers 
Figures 5 7 to 5 9 show the effect of crystal rotation on buoyancy driven con- 
vection m the melt When there is no crystal rotation (case (a) m Figures 5 7 to 
5 9), flow IS completely driven by buoyancy and counter clockwise convective rolls 
are initiated along the crucible wall Superposition of crystal rotation with buoy- 
ancy driven natural convection (case (b) and (c) m Figures 5 7 to 5 9) brings a 
two cell pattern in the flow field As the crystal rotation rate increases, the clock- 
wise convective rolls initiated due to crystal rotation push the natural convective 
currents towards the wall Pumping action of crystal lotation makes the hot fluid 
near bottom move up towards crystal and thus isotherms in the temperature field 

shift upward towards the crystal ^ 

In Figoie 6 10, variation in pull velocity of a continuously growing Nd YAG 

crystal for two diflerent initial melt heights is shown Both the plots show a 
continuous decrease m pull velocity with time until it becomes zero As the 
crystal grows, melt height goes down In these simulations, this lowering in melt 
height IS very-very small as compared to the mcisase m crystal height Thus 
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duiing flow and temperatuie fields in the melt changes almo-it neghgibh v hile 
ladiation and convection losses fiom cpstal stiiface contmnousb increase 
the melt ci3''Stal mteiface melt side heat fin\ leinaiiis almost ‘^amc vhile solid 
side flu\ incieases continuoiisl} with inr leasing cnstal height Tins makes the 
diffeience between the two fluxes deciease continuoush till it becomes zero Since 
late of solidification is piopoitional to this flux unbalance pull \elocitj' decreases 
continuously and finally becomes zeio 

Vaiiation in concentiation field in the melt wutli time foi diffcient cr\stal 
lotations has been showm in Figuies 5 11 to 5 16 Concentiation contours sliowm 
aie foi segregation coefficients, — 0 (Figuies 5 11 to 5 13) and ko = 0 2 (Figures 
5 14 to 5 16) No change m concentration of the melt takes place foi A-o = 1 
Due to scgiegation concentration of the solute inci eases near the melt-civstal 
mteiface and this solute is then ledistnbuted in the melt by con\ection currents 
Concentiation contours show how' solute redistribution is caused b\ complex flow 
patterns in the melt Rate of melt enrichment is maximum for ko = 0 as all the 
solute m the solidifying melt is rejected Rejection foi Aq = 0 2^ corresponds to 
the Nd YAG case 

Figure 5 17 and 5 18 shows ladial distribution of solute just below the melt- 
crystal interface It can be seen that for no crystal rotation, radial variation 
m concentration is maximum with maximum concentration at the center of the 
crystal This is simply a consequence of the purely buoyancy driven flow wffiich 
makes the solute move towards crystal-center Now as the crystal rotation is 
increased concentration profile starts becoming flat due to the centrifugal action 
at the melt-crystal interface 


value of Aq = 0 2 has been 


found m the experiments of Lw [2000] 









Figure 5 7 Temperature (isotherms on the left) and flow field (streamlines on 
the light) in Nd YAG-melt for thiee different crystal rotations (t e , (a) Re = 0, 
(b) Re = 950 and (c) Re = 500) For all thiee ciystal lotations, melt height 
and Grashof numbei are 1 and 1 x 10“^ respectively 
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Figure 5 9 Temperature (isotherms on the left) and flow field (streamlines on 
the right) in Nd YAG-melt for three diffeient crystal lotations {i e , (a) Re = 0, 
(b) Re 250 and (c) Re = 500) Foi all three crystal rotations, melt height 
and Grashof number are 0 5 and 1 x 10^ lespectively 






time, s 


(b) 

Figuie 5 10 Vaiiation m pull velocity of a continuously growing Nd YAG crystal 

foi two diffeient initial melt heights, ? e , (a) 1 and (b) 0 5 In both the cases, 

initial crystal height is 0 1 and Gr = 1 x 10^ Re is Reynolds number for the 
crystal rotation 
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Figure 5 11 Iso-concentiation lines for Nd concent < 
ent values of time, i e , (a) 0 01, (b) 0 05, (c) 0 1, (c 
variable is taken as {C — Co) x 10^ i where Co(— 1) 
Other parameteis aie Gr = 1 x 10'®, Re — 0, segu, 
pull rate = 3 mm/h 


'I dif melt at six diffei- 
'' j ! and (f) 5 Contoui 
il solute concentiation 
' Of fhcient Lq = 0 and 
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Figuie 5 12 (0^0 MdTo t (e)trd (f)1 'cottom 

cut values of time, ' ^ ^ } ^ J initial solute concentration 

variable is taken as (C - Co) x 10 w^ie^ ol ^ ^ 0 and 

Othei paiameters are Gr = 1 x 10 , Ite - mu, sug 
pull rate = 3 uim/h 
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Figuie 5 13 Iso-concentratioa lines for Nd concentration in the melt at six dfc 
ent values of tune, t e , (a) 0 01, (b) 0 05, (c) 0 1, (d) , (e) an . t,on 

vanable is taken as (C - Co) >- 10^ wheie C.(= 1) .s mmal “‘j?; 

Other parameters are Gr = 1 x 10‘, Re = 600, segregation coefficient h - 0 and 

pull late = 3 mm/ll 




Figuie 5 14 I&o-concentration lines foi \M concentiation m the melt at si\ differ- 
ent values of time, z e , (a) 0 01, (b) 0 05, (c) 0 1, (d) 0 5, (e) 1 and (f) 5 Contour 
variable is taken as (C — Co) x 10®, wheie Co(= 1) is initial solute concentiation 
Other paiametois are Gr = 1 x 10'^, Re = 0, segregation coefficient Aq = 0 2 and 
pull rate = 3 mm/h 
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Figure 5 15 Iso~concentration lines for Nd concentration in the melt at si\ diuer- 
ent values of time, i e , (a) 0 01, (b) 0 05, (c) 0 1, (d) 0 5, (e) 1 and (f) 3 Contour 
variable is taken as {C — Cq) x 10®, where C'o(= 1) is initial solute concentration 
Othei parameters aie Gr = 1 x 10“^, Re = 250, segregation coefficient /cq — 0 2 
and pull rate = 3 mm/h 




Iso-concentration lines for Nd concentiatioii in the melt at 
time, J e , (a) 0 01, (b) 0 05, (c) 0 1, (d) 0 5, (e) 1 and (f) 5 
.ken as (<7 - Co) x 10’ where C'o(= 1) is initial solute conc^ 
leters aie Gr = 1 x 10^ Re = 500, segicgation coefficient 
3 = 3 mm /h 
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Figure 5 17 Radial concentration profiles at different times for three different 
values of crystal rotations Other parameters are^ initial melt height = 1, Gr =: 
1 X 10^5 segregation coefficient Lq 0 and pull rate = 3 mm/h 
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Figure 5 18 Radial concentration profiles at diffeient times for three different 
values of crystal rotations Other parameters aie, initial melt height = 1, Gr = 
1 X 10^, segregation coefficient ko = 0 2 and pull late = 3 mm/h 





Chapter 6 

Conclusions and Scope for Future 
Work 

6.1 Conclusion 

The Preconditioned Conjugate Gradient (PCG) algorithm has been tested for 
various benchmark problems m heat transfer Variation in growth rate, melt 
convection and solutal transport have been studied for the growth of Nd YAG 
single crystals by the Czochralski process Following major conclusions have 
been arrived at m the present work 

1 CPU times for the test case of unsteady-nonlinear heat conduction show 
that PCG IS extremelv effective for matrix inversion as compared to Gaus- 
sian elimination and the Gauss-Seidel technique Moreo\er, PCG becomes 
an even more attractive choice on finer grids 

2 PCG is well-suited for matrix inversion m advection-diffusion problems 
The results obtained for Nusselt number in channel flow match the analyt- 
ical solution 

3 Eigenvalue spectrums of the original and preconditioned matrices show that 
preconditioning with incomplete lower-upper (ILU) decomposition is very 
effective Moreover, it is general m the sense that it can be constructed and 
applied for any nonsingular coefficient matrix 

4 In the present simulation for crystal growth, crystal pull rate is varied to 
grow constant diameter crystals The results obtained for puli-velocity -vaii- 
ation show that it is not possible to maintain favorable giowth conditions 
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for a long peuod of time To grov. longer constant diameter crystals, pa- 
lameters such as crucible tempeiature and enclosure temperature need to 
be vaiied 

5 Redistiibution of the segregated solute in the melt tales place because of 
the complex flow patterns in the melt Moreover, it is not possible to giow 
crystals with uniform radial solute distribution without any crystal rotation 

6.2 Scope for Future Work 

The present work has opened up new aieas of exploration Folloving challenges 
are recommended for future researchers 

1 An extensive assessment of various preconditioning techniques for PCG is 
required 

2 Attractive parallel properties of the effecrive PCG algorithm calls for the 
development of appropriate algorithms 

3 An extensive study of the efltect of various parameters on crystal growth is 
needed to provide guidelines for growing high quality single crystals 

4 Exact prediction of solute incorporation into the grown crystals needs a 
microscopic model for solute segregation 
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